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1 Introduction 


This report describes our research on sparse matrix techniques for the Computational Struc- 
tural Mechanics (CSM) Testbed [22] conducted for NASA grant NAG-1-803. Before provid- 
ing a synopsis of the report, we give a brief overview of the work that has been completed 
during the 10-month tenure of the grant. 

A primary objective was to compare the performance of state-of-the-art techniques for 
solving sparse systems with those that are currently available in the CSM Testbed. Thus, 
one of the first tasks was to become familiar with the structure of the Testbed, and to install 
some or all of the SPARSPAK package [2, 13, 14] in the Testbed. 

We began by installing the CSM Testbed on our SUN workstations. We were the first 
site to do this, and it was necessary to collaborate closely with the CSM group at Langley 
in order to resolve some minor problems with the installation procedure. 

A suite of subroutines to extract from the database the relevant structural and numer- 
ical information about the matrix equations has been written. A driver program (proces- 
sor) that employs these routines along with the SPARSPAK library has been written, and 
we have successfully solved all the demonstration problems distributed with the Testbed. 
These codes have been documented, and performance studies comparing the SPARSPAK 
technology to the methods currently in the Testbed have been completed. In addition, 
some preliminary studies have been done comparing some recently developed out-of-core 
techniques with the performance of the Testbed processor INV. 

An outline of the report is as follows. Section 2 contains a brief overview of the CSM 
Testbed software and its usage. This is essentially background material for the uninitiated, 
and can be ignored by those with experience in the usage of the Testbed. 

Since the ultimate goal of sparse matrix research for the Testbed is to enhance the 
performance and capabilities of the Testbed, some knowledge of the methods currently 
employed is essential in the development of better techniques for the Testbed. Section 3 
gives an overview of the sparse matrix techniques currently employed in the CSM Testbed. 
Our presentation is focused on the internal working of the SPAR matrix processors [25]. 

Section 4 describes an interface which we have designed and implemented as a research 
tool for installing and appraising new matrix processors in the CSM Testbed, along with a 
description of a new processor SPK which consists of a subset of SPARSPAK-A [2] and a set 
of subroutines which provide an interface between SPARSPAK-A and the global database 
of the CSM Testbed. A guide for installing the processor SPK in the Testbed is provided 
in Appendix A of this report. The installation dependent modules of this processor are 
listed in Appendix B with comments indicating the changes to be done at a different site. 
A listing of all interface subroutines is provided in Appendix C. 


3 


PRECEDING PAGE BLANK NOT FILMED 


MUNttMtAUX ft m 


BUM ts 


Finally, Section 5 contains results of numerical experiments we performed in solving 
a set of Testbed demonstration problems using the processor SPK and other experimental 
processors. These results are compared with the performance of the SPAR matrix processors 
on the same set of test problems. 

2 The CSM Testbed Software System and Its Usage 

To facilitate our discussion throughout this report, we shall first briefly introduce the con- 
cepts and terminology employed in the Testbed. Since our discussion is conducted primarily 
for the readers who have not used the Testbed before, the readers who are familiar with its 
usage can skip this section. 

The CSM Testbed is a structural analysis system evolving from integrating the SPAR 
finite element code [25] and the NICE data management and command processing utilities 
[4, 5, 6, 7, 26]. The FORTRAN programs for SPAR (Structural Performance Analysis and 
Redesign) were developed in the 1970’s by Lockheed Missiles and Space Company and by 
Engineering Information Systems, Incorporated. The SPAR system uses the finite element 
approach to perform stress, buckling, vibration, and thermal analysis on linear structural 
systems. The NICE (Network of Interactive Computational Elements) system was originally 
developed at Lockheed Palo Alto Research Laboratories to support engineering analyses. 
The major components of the NICE system include a data manager, a command language 
and a command interpreter. Continued effort has been made by the CSM development 
team at NASA Langley and at the Lockheed Palo Alto Research Laboratory to extend 
the analysis capability of the Testbed since the implementation of its initial version (called 
NICE/SPAR). 

The user interface for the Testbed is described in detail in the CSM Testbed User’s 
Guide [24]. The language, directives, interface, global- database manager and input-output 
manager of the CSM Testbed architecture are each documented in references [5, 6, 7, 8, 26]. 
For our purpose we shall simply walk through an example to quickly familiarize the readers 
with the general usage of the Testbed. The example we use is a Testbed demonstration 
problem presented in reference [25]. We shall refer to this example as problem “demol” 
throughout this report. 

The operating environment Our discussion throughout this report refers to the version 
of the Testbed currently operational on a SUN 3/50 workstation running the UNIX 1 
operating system at the University of Tennessee, Knoxville. 

The problem to be solved: The tubular beam shown in Figure 1 is cantilevered at joint 
1 and statically loaded at joint 5. The static solution for a transverse shear load of 
1000.0 and for an axial load of 10000.0 is required. 

1 UNIX is a trademark of AT&T Bell Laboratories. 
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Tube, inner radius = 2.00, outer radius = 2.25 
E = 10. x 10 6 
v = 0.3 

p = 0.101 

a =0.1 x 10 ~ 4 

Figure 1: CSM Testbed Demonstration Problem - Tubular beam. 


User input Edit a file to contain the script in Figure 2. The command stream demon- 
strates how to solve the tubular beam problem in Figure 1 using the NICE command 
language and the SPAR computational modules. 

Comments The problem-oriented Testbed command language is called CLAMP - 
an acronym for Command Language for Applied Mechanics Processors. The 
commands with their leading keyword prefixed by an asterisk are called CLAMP 
directives. They are special commands used to 

- directly access a global database, 

- define command procedures, 

- implement branching and cycling for nonsequential command processing, 

- process macrosymbols in an advanced language construct, 

- request other available services. 

For example, the directive 

*open 1 demol.101 /new 

contained in our script file will create a new library file with the library identifi- 
cation number (LDI) equal to “1” and file name of “demol.101”. 

The SPAR processors are each implemented as a subroutine callable by the 
Testbed executive module. The macroprocessor command to start the execu- 
tion of a processor is [XQT. Therefore, during the execution of the Testbed, the 
command to run the SPAR processor named TAB is 
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[XQT TAB 


The input (user commands and/or data) to a processor are entered after the [XQT 
command according to the requirements of the individual processor. The SPAR 
input syntax and processor requirements are described in detail in reference [24]. 
Since the CLAMP directives may be intermixed with the processor commands in 
the script file, it is worth noting that once the execution of a processor is initiated 
by [XQT, it will begin and continue accepting input until either another [XQT, a 
STOP or a +ST0P is encountered. If a STOP occurs, execution will proceed to com- 
pletion of the processor’s assigned task after which the next command, which can 
be either a CLAMP directive or a macroprocessor command, begins execution. 
A *ST0P terminates execution immediately. Therefore, the user command STOP 
in the sequence 

[XQT SSOL 
STOP 
*T0C 1 

is necessary to ensure that processor SSOL runs to completion before the directive 
♦TOC is processed. 

The modular structure of the Testbed implies that multiple processors are typ- 
ically executed to perform an analysis. These processors communicate through 
a common database consisting of global-access data libraries (GAL) which are 
operated on by the NICE data manager GAL-DBM [26]. Each GAL data library 
may contain multiple nominal datasets. Each dataset is made up of named 
records. The GAL-Processor interface facilities allow the Testbed processors to 
generate, store, locate, and access all of the needed information in the global 
database to perform a required analysis. The table of contents of an active data 
library may be displayed during execution of the Testbed via the CLAMP di- 
rective *T0C. In Figure 3, we display the table of contents for the data library 
“demol.101” (LDL=1) created by executing the script in Figure 2. 

To execute the analysis: Note that on UNIX systems the execution of the Testbed is 
initiated by the first command “time nicespar << \eof” in the script file, where 
“nicespar” is the name of the executable file and we assume that the name of the 
directory where “nicespar” resides has been inserted in the user’s PATH environment 
variable. Note also that “\eof ” is the last entry of the script. Assuming that the name 
of the file containing the script is “demol.com” and that it has been made executable 
with the “chmod” command, the script may be run by typing 

demol.com 
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To print the solutions on an ordinary text files The default output file for the Testbed 
is the standard output on UNIX systems. The command 

demol.com > ft demol.log ft 

thus redirects the output to the log file. The desired static solutions are produced by 
processor SSOL and the actual data are contained in the dataset named STAT .DISP .1.1. 
To print the static solutions on the log file, the SPAR utility processor VPRT may be 
executed after [XQT SSOL. The command to be inserted into the script is 

[XQT VPRT 

TPRINT STAT DISP 1 1 

The output corresponding to this command is displayed in Figure 4. Note that each 
constrained component is flagged with an asterisk by the processor VPRT. 

More details: We shall come back to this example from time to time to provide the details 
which are not needed until our discussion at a later point. 
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time nicespar « \eof 
♦open 1 demol.101 /new 
♦set echo=off 
[xqt TAB 

START 5 

JOINT LOCATIONS 

1 0 0 0 . 

2 0 0 10 . 

3 0 0 20. 

4 0 0 30. 

5 0 0 40. 

MATERIAL CONSTANTS 

1 10.E+6 .3 .101 .IE-4 
BEAM ORIENTATIONS 
11111 . 

E21 SECTION PROPERTIES 
TUBE 1 2. 2. 25 
CONSTRAINT DEFINITION 1 
ZERO 123456 
1 

[xqt ELD 
E21 
1 2 

2 3 

3 4 

4 5 
[xqt E 
[xqt EKS 

[xqt RSEQ 

reset METH0D=1 LJSPRT=1 LADPRT=1 
[xqt TOPO 

reset PRTKMAP=1 PRTAMAP=1 
[xqt K 

reset spdp=2 
[xqt INV 

reset spdp=2 
[xqt AUS 
ALPHA 

CASE TITLES 
1 ’TRANSVERSE LOAD 
2 ’AXIAL LOAD 
SYSVEC 

APPLIED FORCES 

CASE 1: 1=2: J=5: 1000. 

CASE 2: 1=3: J=S: 10000. 

[xqt SSOL 
[xqt GSF 
[xqt PSF 
stop 

♦TOC 1 
\eof 


. Start and time Testbed execution 
. Open data library 
. Do not echo input 

. Macroprocessor command to execute TAB 
. 5 nodes points in beam 
. Direct TAB input 


. Constrain 6 components of joint 1 
. to be zero 

. Define elements 

. Define element connectivity 


. Create element datasets 
. Calculate element intrinsic 
. stiffness matrices 
. Resequence nodes 

. Form maps which guide the assembly 
. and factorization of system matrices 
. Assemble system stiffness matrix 
. Output dataset in double precision 
. Factor system stiffness matrix 
. in double precision 

. Direct AUS input 
. Define load titles for 2 cases 


. Dir-2 load on joint 5 of 1000. 

. Dir-3 load on joint 6 of 10000. 

. Solve for static displacements 
. Compute stresses 
. Print stresses 

. Print Table of contents of library 1 


Figure 2: A runstream for solving problem demol. 
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

+ Library 1 File: demol.101 + 
+ Form: GAL82 File size: 22062 vords lo. of Datasets: 32 + 
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 


•q* 

Date 

Time 

Lk 

Records 

Processor Dataset name 

i 

05:14:88 

17:54:17 

0 

1 

TAB 

JDF1 . BTAB.l . 8 

2 

05:14:88 

17:54:17 

0 

1 

TAB 

JREF.BTAB.2.6 

3* 

05:14:88 

17:54:17 

0 

1 

TAB 

ALTR.BTAB.2.4 

4 

05:14:88 

17:55:17 

0 

3 

TAB 

GMTR . BTAB . 6 . 6 

5 

05:14:88 

17:55:17 

0 

1 

TAB 

ALTR.BTAB.2.4 

6 

05:14:88 

17:55:17 

0 

1 

TAB 

JLOC.BTAB .2.5 

7 

05:14:88 

17:55:17 

0 

1 

TAB 

MATC.BTAB.2.2 

8 

05:14:88 

17:55:17 

0 

1 

TAB 

MREF.BTAB .2.7 

9 

05:14:88 

17:55:17 

0 

1 

TAB 

BA. BTAB. 2. 9 

10 

05:14:88 

17:56:17 

0 

1 

TAB 

C0B..1 

11 

05:14:88 

17:56:17 

0 

1 

TAB 

QJJT. BTAB.2. 19 

12 

05:14:88 

17:56:17 

0 

1 

ELD 

DEF.E21.1 .2 

13 

05:14:88 

17:56:17 

0 

1 

ELD 

GD.E21.1.2 

14 

05:14:88 

17:56:17 

0 

1 

ELD 

GTIT.E21 .1.2 

IB 

05:14:88 

17:56:17 

0 

1 

ELD 

DIR.E21 .1.2 

16 

05:14:88 

17:56:17 

0 

1 

ELD 

ELTS.BAME 

17 

05:14:88 

17:56:17 

0 

1 

ELD 

ELTS.BBOD 

18 

05:14:88 

17:56:17 

0 

1 

ELD 

ELTS.ISCT 

19 

05:14:88 

17:56:17 

0 

1 

ELD 

IS 

20 

05:14:88 

17:56:17 

0 

4 

E 

E21 . EFIL.l . 2 

21 

05:14:88 

17:56:17 

0 

1 

E 

DEM.DIAG 

22 

05:14:88 

17:56:17 

0 

1 

RSEQ 

JSEQ.BTAB.2.17 

23 

05:14:88 

17:56:17 

0 

1 

TOPO 

KMAP. .9.3 

24 

05:14:88 

17:56:17 

0 

1 

TOPO 

AMAP. .9.3 

25 

05:14:88 

17:56:17 

0 

1 

K 

K. SPAR. 36 

26 

05:14:88 

17:56:17 

0 

5 

IBV 

MV.K.l 

27 

05:14:88 

17:56:17 

0 

2 

AUS 

CASE. TITL. 1.1 

28 

05:14:88 

17:56:17 

0 

2 

AUS 

APPL.FORC.l .1 

29 

05:14:88 

17:56:17 

0 

2 

SSOL 

STAT.DISP.1.1 

30 

05:14:88 

17:56:17 

0 

2 

SSOL 

STAT.REAC.1.1 

31 

05:14:88 

17:56:17 

0 

4 

GSF 

STRS.E21 .1.1 

32 

05:14:88 

17:56:17 

0 

4 

GSF 

STRS.E21.1.2 


++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 


Figure 3: Table of Contents of Library 1. 
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** BEGIN VPRT ** DATA SPACE= 600000 WORDS 
1STATIC DISPLACEMENTS. ID= 1/ 1/ 1 

TRANSVERSE LOAD 


OJOIHT 

1 

2 

3 

4 

5 

6 

1 

0. OOOe+OO* 

0. OOOe+OO* 

0. OOOe+OO* 

0. OOOe+OO* 

0. OOOe+OO* 

0. OOOe+OO* 

2 

0 . OOOe+OO 

0.2506-01 

0. OOOe+OO 

-0.463e-02 

0. OOOe+OO 

0. OOOe+OO 

3 

0. OOOe+OO 

0.8976-01 

0. OOOe+OO 

-0. 793e-02 

0 . OOOe+OO 

0. OOOe+OO 

4 

0. OOOe+OO 

0.1816+00 

0. OOOe+OO 

-0.992e-02 

0. OOOe+OO 

0. OOOe+OO 

5 0. OOOe+OO 0 . 285e+00 

1 STATIC DISPLACEMENTS. 

AXIAL LOAD 

0. OOOe+OO 

-0. 106e-01 

0. OOOe+OO 
ID= 1/ 

0. OOOe+OO 
1/ 2 

OJOINT 

1 

2 

3 

4 

5 

6 

1 

0.0006+00+ 

0. OOOe+OO* 

0 . OOOe+OO* 

0. OOOe+OO* 

0. OOOe+OO* 

0. OOOe+OO* 

2 

0.0006+00 

0. OOOe+OO 

0.300e-02 

0. OOOe+OO 

0. OOOe+OO 

0. OOOe+OO 

3 

0.0006+00 

0. OOOe+OO 

0.599e-02 

0. OOOe+OO 

0. OOOe+OO 

0. OOOe+OO 

4 

0.0006+00 

0. OOOe+OO 

0 . 899e-02 

0. OOOe+OO 

0. OOOe+OO 

0. OOOe+OO 

S 

0.0006+00 

0. OOOe+OO 

0.1206-01 

0. OOOe+OO 

0. OOOe+OO 

0. OOOe+OO 


EXIT VPRT CPUTIME= O.S I/0(DIR,BUF)= 0 0 


Figure 4: The contents of dataset STAT.DISP.1.1. 
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3 The CSM Testbed Matrix Processors 


Reference [23] contains a set of logic flowcharts developed for the key subroutines of each 
of the SPAR matrix processors TOPO, K, INV, SSOL and AUS. These charts together with 
the commented FORTRAN source code are very helpful in our understanding of the sparse 
matrix techniques currently employed in the Testbed. In this section, we shall attempt to 
describe the algorithms and data structures which are implemented by the processors INV 
and SSOL. 

3.1 The Basic Algorithms 

The factorization algorithm: Processor INV applies a specialized Gaussian elimination 
scheme to factor a sparse symmetric matrix K into LDL T , where A is a unit lower 
triangular matrix and D is a diagonal matrix. This algorithm is numerically stable if 
the matrix K is also positive definite, which is the case when K is the system stiffness 
matrix. The basic algorithm can be easily described for a dense symmetric matrix A 
as follows. We assume that A is of dimension n x n. Let us denote the elements of A 
and M - L T as a,j and m tJ -, where 1 < i < n and i < j < n, and D = {di , d 2 , . . . , d n }. 
Note that each off-diagonal a tJ - is overwritten by rriij and that each an is overwritten 
by d^ 1 if the algorithm presented in Figure 5 is successfully executed. Algorithm I 
assumes that the a tJ - elements are stored row by row. 


Algorithm I The basic LDL T factorization scheme 

for i *— 1, 2, . . . , n do 
if an = 0 then 
quit 
else 

a u * 1 Jan 

for k <— i -f 1, . . . , n do 
for j <— n do 

flfcj * ~ 1 &kj ^ ^ Q* ij 

for k <— t + 1, . . . , n do 

Oj k < 


Figure 5: Computing D~ l and M — L T factors of A. 
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The following features of the algorithm above will be exploited in its sparse imple- 
mentation. 

1. To compute D" 1 and the off-diagonal elements of M — L T , the elements stored 
and accessed are those on the diagonal and in the upper triangular part of A. For 
example, when n = 5, the algorithm performs the transformation in Figure 6. 


a l,l a l,2 

a l,3 

Ol,4 

a l,5 


d.- 1 m i ,2 m li3 m 14 m li5 

a 2,2 

02,3 

a 2,4 

a 2,5 


dj 1 m 2,3 m 2 , 4 m 2 ,5 


«3,3 

03,4 

03,5 


dj 1 " 13,4 m 3 5 



04,4 

04,5 


d- 1 m 4,s 




05,5 


ds 1 


Figure 6: Overwriting AbyD~ l and M = L T . 

2. The dij s which have been overwritten by the elements of D~ l and M — L T 
will not be needed in the remaining elimination stages. In particular, during the 
i th elimination stage, the elements accessed and modified are confined to row i 
through row n as shown in Figure 7 for i = 3 and n = 5, where 0 represents 
elements which are not accessed. 




® 

® 


® ® 

® 

0 0 


® 

® 

® 


® 

® 

0 0 


03,3 

03,4 

03,5 



dj 1 

m 3,4 m 3,5 



04,4 

04,5 




^ 4,4 ^45 




05,5 




^ 5,5 


Figure 7: LDL T factorization of A - the third stage. 
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Solving the triangular systems: Since Algorithm I stores the factors D -1 and M = L T , 
we shall describe the solution scheme in terms of these two factors. Both of the forward 
and backward substitution schemes presented below access the elements of the factor 
M row by row. 

Step 1. Forward substitution scheme (Solve M T y = b), 

for i <— l,...,n do 
Vi *- k 

for i = i+l,...,ndo 
bk — bk — m ik * yi 

Step 2. Backward substitution scheme (Solve Mx = D^y). 
for p <— 1 , . ..,n do 

i <r~ 71 — p + 1 

s * — 0 

for j = i + 1, . . .n do 
s <— s + rriij * Xj 
X{ 4- dr 1 *y { - s 


3.2 The INV Implementation 

In this section, we shall discuss in various degrees of details the following aspects of the 
sparse factorization scheme implemented by the processor INV. 

1. The algorithm - a block LDL T factorization scheme. 

2. Memory requirement. 

3. Data structures. 

4. The handling of zero constraints. 

5. The handling of nonzero constraints. 

6. Data archived to the global database. 

A block LDL t factorization scheme: The processor INV has tailored Algorithm I to 
perform an out-of-core block LDL T factorization of large sparse matrices arising in 
the finite element analysis of structural mechanics problems. Before we describe the 
INV implementation of this scheme, let us first explain the block LDL T algorithm 
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by applying it to a dense symmetric matrix A in block form. To be specific, let us 
consider the 2x2 block matrix in Figure 8, where Ai t i = {a£J}, Ai t2 = 

A 2 ,2 = with = a^, aff = and 1 < k,j < 3. 


A = 




4:1 

a (‘) 

a i,2 

a (>) 

° i,j 

42 

a (”) 

° 1,2 

42 \ 



4?1 

a (») 

a 2,2 

fl (0 

“ 2,3 

4!? 

a ( ”> 

a 2,2 

4"’ 

> 11,1 

Al ,2 \ 

a «> 

° 3,1 

a (0 

° 3,2 

4‘J 

42 

a <«> 

“ 3,2 

a <«> 

“ 3,3 

a t 

> 12,2 ) 

“l,l 

a (*') 

a 2,l 

a (-o 

“ 3,1 

i “i:: 1 

W <T 

*(•») 

“ l,3 



42 

a (*0 

a 2,2 

«(«) 

“ 3,2 

i “4’ 

4‘r) 

“ 2,2 

4!;’ 



\ <1 

«(*) 

a 2,3 

a (0 

a 3,3 

i 4? 

“ 3,2 

42 / 


Figure 8: Partitioning symmetric A into four 3x3 blocks. 

The block LDL T scheme works in the following manner. 

Step 1. Apply Algorithm I to matrix Ai t i to perform the following transformation. 


/ flW ag ag \ f 1/4° 


a (’) a (0 
a 2,2 u 2,3 


x 3,3 


4,2 m l‘,3 \ 


l/d 


(l) m 2,3 


i/4° / 


In other words, at the end of step 1, we have in fact zeroed out the nonzeros 
in the lower triangular part of A lt \ and stored the multipliers Uj = mji in the 
upper triangular part of 
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Step 2 . Apply the multiplers to the Ai t2 block as if LU ( U = DL T ) decompo- 
sition were applied to reduce (Ai t h Aij) to an upper trapezoidal matrix. That 
is, the Ai t 2 block is overwritten by the resulting {uj^} of the following transfor- 
mation. 



Step 3. Zero out the block A [ 2 implicitly by applying the multipliers directory to 
block A 2 > 2- The multipliers can be computed on the fly from 



The A 2t2 block is then updated to be A 2i2 = {a£j^} which is obtained by the 
following computation. 



Since the A 2>2 diagonal block is symmetric, only the upper triangular part of 
A 2,2 is updated in the actual computation. 

Step 4. tijjjj «“ u kj/ d k\ for Vfc > 3- 
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Note that the transformations accomplished by the above four steps can be expressed 
with respect to the block upper triangular part of the given matrix as follows. 


/ <* 1,1 

a l ,2 

<* 1,3 

a l ,4 

a l ,5 

<* 1,0 \ 


dr 1 

*** 1,2 

*** 1,3 

7711,4 

*** l >5 

*** 1,6 

\ 

<* 2,1 

<* 2,2 

<* 2,3 

<* 2,4 

<* 2,5 

<* 2,6 


<* 2,1 


*** 2,3 

*** 2,4 

*** 2,5 

*** 2,0 
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The final step: For this particular example, the factorization is completed after 
transforming 


(a 4)4 a 4i 5 a 4i6 

65,5 «5,6 

a <>,6 


1 d 4 x m 45 m 4 , 6 \ 

ds 1 rn s , e 
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by Algorithm I. 

The output matrix: The coefficient matrices of the resulting triangular systems, 
namely M T y — b and Mx = D~ x y , are available from the output matrix given 

by 


l d r 1 


m 1,2 
d; x 


***1,3 

m i,4 

***1,5 

m i,6 

***2,3 

m 2,4 

***2,5 

m 2 ,6 

dz 1 

m 3i 4 

***3,5 

m 3,6 


V 

***4,5 

m 4i6 




m 5,6 




dr 1 
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The following observations may be made on the block LDL T factorization scheme 
described above. 


1. The elements in the lower triangular part of the diagonal Aij blocks are not 
accessed during the process of computing the D _1 and L T = M factors. 

2. The block of rows which have been updated to contain the d j^’s and m^j’s of 
the factors are no longer needed in the future stages of elimination. 

3. Observe that the updating of A 2l 2 block in Step 3 can be reformulated as follows. 
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and 


s(**) 

° 1,1 

a (iv) 

°1,2 

*(»») \ 
°1,3 
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Therefore, if the elements of A 2| 2 are not available in memory at the time the first 
block row (A^i, Ai | 2 ) is being processed, the modifications can be accumulated 

in the which are later added to the respective elements of A 2|2 when 

they are read into memory. 


Although it appears straightforward to generalize the block LDL T scheme to a sym- 
metric sparse block matrix such as the example given in Figure 9, where each is 
a dense square matrix of some uniform dimension, an efficient implementation of the 
sparse block LDL T scheme requires sophisticated data structures. 


k 1a 

^1,2 



K llS 
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*2,« 



^3,4 

*3,5 

*3,« 



*4,4 


K 4,8 




*5,5 

*5,6 







Figure 9: Upper triangular block structure of a symmetric sparse matrix K. 


Memory requirement: Suppose that the matrix K in Figure 9 is stored out-of-core and 
the rows of K are to be read into memory one block row (i.e., JDF rows if JDF is 
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the dimension of each submatrix) at a time. In order to factor the first block row 
K\faK\$) and store the modifications to be applied to the blocks #2,2, ^2,5 
and K $ } 5 later, we need memory space to store the blocks in Figure 10 as well as 
the indexing overhead incurred by the data structures employed. In order to proceed 


*1,1 

*1,2 

*1,5 


-$2,2 

-^2,5 



*^5,5 


Figure 10 : The storage needed for processing Kl.s) 

with the factorization of the second block row (If 2,2, 1^2,3, If 2|5 , 1^2,6), enough working 
space must be available to accommodate the blocks in Figure 11 . To minimize the 
memory requirement, processor INV actually re-uses the space occupied by blocks K 11, 
K 12 and ff 1(5 to accommodate the blocks needed for the current elimination stage, 
assuming that the factors of (Ifij, Ifi ( 2, ^1,5) have been archived to the database. 
The block submatrices needed to remain in memory for each of the next four stages 


^2,2 

*3,3 

-^2,5 

*2,8 


S3, 3 

S3 ,5 

S3 ,6 



S5 ,6 

S5 ,5 




S 6| 6 


^ 2,2 

^ 2,5 


* 2,2 + S 3 ,2 

*2,5 + S 2t 5 


Figure 11 : The storage needed for processing (K 2)2 , ^2,3, #2,5, #2,6) 

are depicted in Figure 12 . Observe that although if 45 block is null in /C, it is to be 
filled in the third elimination stage. Therefore, the space to accommodate S 4)5 block 
must be allocated. Fortunately, the fill-in locations can be determined prior to the 
numerical factorization phase. With the fill-in information available, the maximum 
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Figure 12: The storage needed to process #3,*, K 4,*, ifs,* and Ke,** 

number of submatrices ever needed to be in memory can also be determined. As 
far as the indexing overhead is concerned, a simple and effective strategy is to store 
one pointer for each submatrix assuming that the elements within each submatrix are 
stored in consecutive locations. Using this indexing strategy, the number of pointers 
required to be in memory for each particular elimination stage is equal to the number 
of submatrices to be present. 

Data structures: The data structures employed by the Testbed matrix processors can 
again be more easily explained using our block 2x2 example given in Figure 8. 

Data structure of the input coefficient matrix: Processor INV assumes that the 
block upper triangular part of the coefficient matrix 



is stored out-of-core in a block-row-oriented manner. That is, the data of the 
blocks are stored in a one dimensional array following the block sequence as 
depicted in Figure 13. 



Figure 13: The block sequence of input matrix. 

Within each Aij block, the elements are stored column by column. For example, 
the elements of the A\ t 1 block are stored following the sequence in Figure 14. 
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<*l,l 

<*2,1 

<*3,1 

<*1,2 

<*2,2 

<*3,2 

<*1,3 

<*2,3 

“3,3 


Figure 14: The element sequence of block -4i,i. 

Data structure of the working array St Except for the first block row of the in- 
put matrix, the data retrieved from the buffer for each block row are necessarily 
updated by adding the modifications {s,j} accumulated in a working array S. 
Therefore, it is not surprising that the elements within each JDFxJDF ( JDF= 3 
in our example) Sij block are stored in the same manner as the input Aij block. 
For our example, the diagonal block A 2 ,2 must be updated by £ 2,2 before it can 
be factored. Suppose that the respective addresses of these two blocks in the 
buffer and working array are given by the pointers KMAP(IX) and AMAP(JX) 
as depicted in Figure 15. 



AMAP(JX) 


Figure 15: Indexing the buffer and working arrays. 
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The integration of ^ 2,2 into the working array is accomplished by the following 
segment of FORTRAN statements. 


c 

C GET POINTER TO THE BUFFER 

C 

LKSUB * KMAP(IX) 

C 

C GET POINTER TO WORKSPACE STORAGE FROM AMAP 

c 

K = AMAP (JX) 

DO 100 J * 1, JDF 
DO 200 I = 1, JDF 

S (I , J ,K) « S(I,J,K) + A (LKSUB) 

LKSUB = LKSUB + 1 
200 CONTINUE 

100 CONTINUE 


For a general sparse matrix, because the data stored in the working array S are 
dynamically changed by accommodating new data in the space occupied by data 
which have been written out to the database, the Sij blocks corresponding to the 
consecutive Aij blocks in the buffer array may not be neighbors in the working 
storage. To integrate NSUBS ( NSUBS > 1) Aij ’ s into 5, the starting address 
of each Sij must be computed from AMAP each time, resulting in the revised 
code segment. 


c 

c 

c 


c 

c 

c 

c 


GET POINTER TO THE BUFFER 


LKSUB = KMAP(H) 

DO 300 ISUB = 1, NSUBS 


GET POINTER TO WORKSPACE STORAGE FROM AMAP 
JGAP IS KNOWN FROM THE DATA FORMAT OF AMAP 


JI * JX + JGAP 
K = AMAP(JX) 

DO 100 J » 1, JDF 
DO 200 1*1, JDF 

S(I,J,K) = S(I, J ,K) + A (LKSUB) 
LKSUB = LKSUB + 1 
200 CONTINUE 

100 CONTINUE 

300 CONTINUE 
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Conversion of the input data structure: Note that the data structure described 
above for the input buffer and working array is in fact the output format of 
the processor which assembles the system stiffness matrix from the finite ele- 
ment model. Since the block LDL T factorization scheme and the following for- 
ward/backward substitution algorithms are row-oriented, the properly updated 
JDFxJDF submatrices of the current block row are copied from S into another 
one- dimensional array B, where the data are stored row by row with respect to 
the global matrix. For example, assuming that the dimensions of 5 and B are 
declared as S( JDF, JDF,*) and B(J DF, CONRNG, *) , the following FORTRAN 
statements retrieve the first row of (Ai t i, Ai^), i.e. {a X|1 , a 1>2 , • • • , ai^}, from S 
and store them in the consecutive locations in B . 


K INDEXES THE CURRENT ROW IN B 


K « 1 


M INDEXES THE CURRENT ROW IN S 


M = 1 


OBTAIN THE NUMBER OF BLOCKS IN CURRENT ROW 


CONRNG = 2 

DO 100 J= 1, CONRNG 


ASSUME THAT THE LOCATION OF THE CURRENT 
BLOCK IN S CAN BE OBTAINED FROM SUBMAP (J) 


LKSUB = SUBMAP (J) 

DO 200 I = 1, JDF 

B (I , J ,K) = S(M,I,LKSUB) 

CONTINUE 
CONTINUE 

Since the modifications computed from B are to be accumulated into S for up- 
dating the input matrix in the future stages of the elimination process, the di- 
mensioning of B as B( JDF, CONRNG,*) in parallel with the dimensioning of 
S is desirable. The conversion of index from S to B, or vice versa, for each 
element can thus be easily expressed in FORTRAN as demonstrated in the 
above code segment. However, there are other times the code would be much 
cleaner by viewing B as a two dimensional array declared as B(JDFCON,*), 
where JDF CON -JDFxCON . The technique which the processor INV uses to 
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index the same array in either way is to declare two formal parameters, namely 
B(JDF,CONRNG,*) and BB(JDFCON f *) in the subroutine which does the fac- 
torization, whereas the actual parameters corresponding to B and BB in the 
calling sequence are identical . With this trick, B and BB in the subroutine refer 
to the same actual parameter and the programmer can work with either B or 
BB according to his need to access the data in a particular pattern. 

Handling zero constraints: Processor INV handles zero constraints by ignoring the cor- 
responding rows in the process of transferring data from 5 to B . That is, if the 
unknown X{ = 0, then row i will not be copied to B. For example, if it is known that 

x 2 = 0, then only row 1 and row 3 in (Ax,x, j4x, 2) would be copied to array B. The 

actual transformation of A\ i2 ) is carried out in B as shown below. 

( ai t i a 1>2 ai j3 ax, 4 ax,s ax,e \ > f d± l m \ y2 *ni ,3 nix, 4 ui 1>5 m lj6 

a 3f x a 3j2 a 3 , 3 a 3> 4 a 3)5 a 3>6 ) \ a 3> i a 3)2 d^ 1 m 3j 4 7713,5 m 3|6 

Consequently, row 3 in S becomes row 2 in J 5, i.e., it is possible that K<M in our 
sample code segment. 

Handling nonzero constraints: Processor INV handles nonzero constraints by ignoring 
the corresponding rows in the factorization process. For example, if it is known that 
z 3 = ^ 0 in addition to x 2 = 0, the transformation of (Ax,i, ^, 2 ) by processor INV 

will not affect row 3, i.e., 


( «l,l a l,2 <* 1,3 a l ,4 <*1,5 « 1,6 

\ a 3 ,x <*3,2 a 3 , 3 <*3,4 a 3j5 a 3f e 


d± 1 Ulx,2 Ulx ,3 m l ,4 *Ul ,5 m l,6 

a 3 ,i <* 3 , 2 a 3}3 a 3) 4 a 3) 5 a 3>6 


Elements archived: Write out to database those elements of BB which are needed for 
the subsequent use by processor SSOL in effecting the forward/backward substitution 
process. For example, assuming x 2 = 0, and x$ = u 3 ^ 0, the output elements 
resulting from factoring the (Ax,i, ^ 1 , 2 ) block are given by 

( d ^ 1 m lt2 m 1)3 mx |4 m 1>5 m 1( e 

V a 3,3 a 3j4 a 3|5 a 3(6 


3.3 The SSOL Implementation 

Input Data: Processor SSOL retrieves from the database the factors archived by processor 
INV. For our example of the block 2x2 matrix, assuming that the constraints are 
x 2 = 0 and x 2 = u 3 ^ 0, the data given below are stored in a row-oriented manner in 
the database. 

df m 1|2 mx, 3 uix,4 rnx,5 mi ,0 ^ 

<*3,3 a 3 , 4 <*3,5 <*3,0 

dl 1 m 4 , 5 m 4) e 

ds 1 m 5i Q 
d, s* 1 
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In addition to the factors, the right-hand side vector / and the nonzero- constraint 
vector u are also available in the database. 


Handling constraints: In essence, processor SSOL simply adapts the forward/backward 
substitution schemes we presented for Algorithm I to solve the following triangular 
systems, which are to be implicitly formed from the data retrieved. 


/ 1 
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d; 1 
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1 7^5,6 


*4 


d; 1 
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1 ) 


*5 


V d^ J 
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In particular, the SSOL implementation takes advantage of the following observations. 


1. The equations corresponding to zero constraints can be ignored in the forward 
substitution phases. 

2. The coefficients of the equations corresponding to nonzero constraints are needed 
to adjust the right-hand side vector in the forward substitution phase. 

3. If the solution vectors contain the constraints, the equations corresponding to 
constraints (either zero or nonzero) can be skipped in the backward substitution 
phase. 


Output Solutions: The computed aq-j’s are written out to the global database. 


3.4 Other Relevant SPAR Processors 

In order to briefly introduce the functions of other relevant SPAR processors, and give the 
readers some idea how they may be used to perform an analysis, we found that the following 
information available in The CSM Testbed User’s Manual [24] useful. Given below is a list of 
processors together with comments on their individual functions. In addition, the ordering 
of the processors in the list serves as a template for performing the linear static analysis, 
which is one of the simplest types of analysis which can be performed with the Testbed. 

1. Processor TAB. Define joint locations, constraints, reference frames, and possibly ma- 
terial and section properties. Material and section properties may be defined using 
either processor TAB or processors AUS and LAU (Steps 2 and 3). 
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2. Processor AUS. Build tables of material and section properties if the facilities in pro- 
cessor TAB were not used. 

3. Processor LAU. Form constitutive matrix if material and section properties were not 
input in processor TAB. 

4. Processor ELD. Define elements. Element definitions include element connectivity, 
element material reference frame number, element material and section type numbers. 

5. Processor E. Initialize element datasets; create the dataset which will contain all im- 
portant element information (e.g., intrinsic coordinates, element-to-global transfor- 
mations, intrinsic stiffness matrices). 

6. Processor EKS. Calculate element intrinsic stiffness matrices. 

7. Processor RSEQ. Resequence nodes for minimum total execution time. 

8. Processor TOPO. Form maps which guide the assembly and factorization of system 
matrices. 

9. Processor K. Assemble global (system) stiffness matrix. 

10. Processor INV. Factor system stiffness matrix. 

11. Processor AUS and EQNF. Create applied nodal loading. If element loading is applied, 
Processor EQNF must be executed to calculate equivalent nodal loading. 

12. Processor SSOL. Solve for static displacements. 

13. Processor GSF. Calculate element stress resultants. 

14. Post-process using any of the following processors: VPRT, PSF, PLTA, PLTB, PLOT, 
CONT, T2PT. 
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4 


Developing New Matrix Factorization Processors 

4.1 General Considerations 

We have described in detail in Section 3 the internal working of processors INV and SSOL. 
The former performs the out-of-core LDL T factorization of a sparse matrix in block form, 
and the latter solves the resulting triangular systems by forward and backward substitution 
schemes. The following considerations have prompted us to investigate alternative sparse 
factorization schemes. 

1. The techniques employed by INV are particularly tailored to the large sparse linear 
systems arising in the structural models. The models considered are composed of 
finite elements connected at specified joints. Each joint can have three translational 
and three rotational components of deflection, totaling a maximum of six degrees of 
freedom per joint. The system stiffness matrix is stored and operated on as an array 
of JDFxJDF submatrices, where 3 <JDF< 6 is the maximum number of degrees of 
freedom per joint in the model of a particular problem. However, in general the joints 
need not have the same number of degrees of freedom. This storage scheme thus 
necessitates storing dummy data - an identically zero row for each missing degree 
of freedom at each joint. Although the factorization scheme only operates on the 
non-null submatrices and some operations on the dummy rows are skipped by the 
processor INV, it does not fully exploit the sparsity within each submatrix. While 
this strategy is understandably very efficient if uniform degrees of freedom per joint 
prevail, it may not best suit the models with drastically varied degrees of freedom, 
which is not uncommon in finite element modeling applied to disciplines other than 
mechanical structural analysis. 

2. As described in Section 3, the data structures employed incur the index overhead of 
one pointer per submatrix for all submatrices occurring in each elimination stages. 
Therefore, the index overhead is proportional to the number of submatrices instead of 
the size of them. Consequently, while the primary storage for the system stiffness ma- 
trix and the factors is reduced for models with fewer degrees of freedom, the secondary 
storage for their indices may remain the same and could become a significant part of 
the total storage. Furthermore, unlike the working storage which is determined by 
the maximum number of submatrices which ever occur during the entire factorization 
process, the addresses of the submatrices are repeatedly stored for each elimination 
stage. 

3. The system stiffness matrix, the factors and their respective indexing information are 
each stored in separate datasets in the global database. The datasets are read into 
memory or written out to the database one record at a time. The choice of record 
length determines the number of disk read/write operations and the required buffer 


26 


space. While the maximum record length of a dataset is restricted by the available 
buffers, the minimum record length must be long enough to contain all of the items 
which are needed to completely process one entire row of submatrices. Therefore, 
the processor INV can perform in-core factorization if each record of each dataset 
contains all information needed to complete the entire factorization process. In that 
case, the in-core storage is required to accommodate at least one copy of the system 
stiffness matrix, one copy of the factors along with the indexing information needed 
for all elimination stages, and a working array of the same size as needed in the 
out-of-core case. Since some other out-of-core sparse factorization schemes currently 
available perform in-place factorization and are readily adapted to performing in-core 
factorization, it appears worthwhile to compare their performance in both in-core and 
out-of-core cases. 

4. When applying the out-of-core block LDL T scheme as implemented by the processor 
INV to a dense matrix, its advantage of reducing memory requirement disappears 
because the working array for the first elimination stage must contain the entire upper 
triangular part of the stiffness system matrix. 

5. The possible ill- conditioning of the system stiffness matrix is not detected by the 
current Testbed software. 

4.2 The Design of an Interface 

It is apparent from our earlier discussions that the format of the datasets is directly con- 
nected to the factorization scheme currently employed in the Testbed. It is thus likely that 
the particular arrangement of data items in the datasets may not be compatible with the 
data-accessing pattern of the other factorization algorithms to be considered. In order to 
evaluate the performance of alternative sparse factorization schemes in the Testbed without 
redesigning the database at a time when the scheme of choice is not certain yet, we have 
devised a set of subroutines which serve as an interface between the global database of the 
Testbed and SPARSPAK-A [2]. Although some components of the interface are specific for 
SPARSPAK-A, we hope that its overall design and the availability of some utility modules 
will prove to be useful in adapting the interface to work with other sparse matrix solvers. 
A few words about the capabilities of SPARSPAK-A are in order. 

4.2.1 SPARSPAK-A: Waterloo sparse linear equations package 

In this section we briefly review the important features of SPARSPAK-A, which is a package 
of Fortran programs designed to efficiently solve large sparse systems of linear equations by 
direct methods. The structure and use of the package are described in the SPARSPAK-A 
User’s Guide [2]. The collection of algorithms implemented by SPARSPAK-A and their 
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storage schemes are discussed in reference [12]. Although we shall consider only symmetric 
positive definite problems here, the actual package handles both symmetric and unsymmet- 
ric problems subject to the condition that the matrix structure is symmetric and that row 
and/or column interchanges are not required to maintain numerical stability. To solve a 
sparse symmetric positive definite linear system 

Ax — b , 

the user and SPARSPAK-A interact through the following steps: 

Step 1 . The user supplies the nonzero structure of A to the package using a set of subrou- 
tines described in Section 2.2 of reference [2]. 

Step 2. The package finds a “good” ordering (permutation P) for A, and allocates stor- 
age for the triangular factorization of PAP T = LL T , as described in Section 2.3 of 
reference [2]. 

Step 3. The user supplies the numerical values for the matrix A to the package, as described 
in Section 2.4 of reference [2]. 

Step 4. The package factors PAP T into LL T , as described in Section 2.5 of reference [2]. 

Step 5. The user supplies numerical values for 6, as described in Section 2.4 of reference [2]. 
(This step may come before Step 4, and may be intermixed with Step 3.) 

Step 6. The package computes the solution by solving Ly = Pb and L T z = t/, and then 
setting x = P r z, as described in Section 2.5 of reference [2]. 

Step 7. The user may call a subroutine to obtain an estimate of the relative error in x 
as well as the inverse of the condition number of A if so desired. The subroutine is 
described in Section 2.6 of reference [2]. 

The names of the subroutines available for reordering a symmetric matrix in Step 2, 
together with the algorithms they implemented, are listed in Table 1. Corresponding to 
each ordering choice in Step 2, a different set of subroutines are provided for Steps 3, 4, 6 
and 7. The subroutines used in Steps 1 and 5 are, however, independent of the ordering 
methods. 

In the context of comparing the performance of the SPARSPAK-A factorization algo- 
rithm with that of the Testbed processor INV, we should note the following. Firstly, the 
coefficient matrix A will have been ordered differently because the ordering algorithm in 
the Testbed is applied to the joints in the finite element model before the system stiffness 
matrix is assembled, whereas SPARSPAK-A reorders the coefficient matrix itself. Since 
associated with each joint in the finite element model is a dense JDFxJDF submatrix, 
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Table 1: SPARSPAK-A ordering choices. 


SPARSPAK-A 

Subroutine 

Ordering algorithm 

ORDRA1 

Reverse Cuthill-McKee ordering [21] 

0RDRA3 

One-way Dissection ordering [9] 

0RDRB3 

Refined quotient tree ordering [10] 

0RDRA5 

Nested Dissection ordering [11] 

ORDRB5 

Minimum Degree ordering [17] 


the resequencing of the joints relocates the submatrices (as a whole) in the system stiffness 
matrix. On the other hand, since the ordering algorithms in SPARSPAK-A are applied to 
the structure of the assembled system stiffness matrix, the zeros within each submatrices 
(due to constrained variables or dummy rows) may be exploited and the resulting matrix 
may not be in block form. 

Secondly, the Cholesky factorization scheme and the upper/lower triangular system 
solvers implemented by SPARSPAK-A do not handle constraints or dummy rows (rows of 
zeros). It is therefore necessary to adjust both the system stiffness matrix and the right-hand 
side before the nonzero structure and the numerical values are input to SPARSPAK-A. In 
the current version of Testbed, while the constraint information is available in a designated 
dataset, the dummy rows can be detected only by reading the assembled system stiffness 
matrix. The implication is that the system stiffness matrix has to be examined twice - 
once for determining its “adjusted” nonzero structure (needed in Step 1), and once for 
retrieving its numerical coefficients (needed in Step 3). We consider the way we handle 
the dummy rows as an interim measure until the dataset format of the generic element 
processor is available. It is expected that the generic element processor will neither assume 
uniform degrees of freedom nor store dummy data. Complete details on adjusting the 
nonzero structure and the numerical values for input to SPARSPAK-A are given later in 
this section. 

Thirdly, SPARSPAK-A employs a particular version of the Cholesky factorization algo- 
rithm. Since this version of the algorithm computes the Cholesky factor one column at a 
time and the part of the matrix remaining to be factored is not accessed during the scheme, 
it is commonly referred to as the “Column- Cholesky” algorithm. Depending on how the 
modifications to each designated column are accumulated, the Column- Cholesky algorithm 
can be described in two different forms. Given in Figure 16 is the commonly known scalar- 
product form. These formulas can be derived directly by equating the elements of A to the 
corresponding elements of the product LL T . 
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for j 
l 


1,2, ...,n do 


3,3 \j a 3y3 

for i <- j + 1J + 2,...,n 


4,j ( a *,j ftjj 


Figure 16: The scalar-product Column- Chole sky Factorization Algorithm. 

The vector-sum Column- Cholesky algorithm described in Figure 17 is an alternative 
formulation which avoids explicitly forming the individual inner products. 


for j = 1,2, ...,n do 

for k — 1, 2 , . . . , j — 1 do 


( a jj 

y 

( a n > 

-bk 

hk \ 

\ a nj 

) 

\ a n j ) 


V Lk J 


hi \/ a Jj 

for fc = j + 1, j + 2, . . . , n do 

h j <*“ a kj/ a jj 


Figure 17: The vector- sum Column- Cholesky Factorization Algorithm. 

SPARSPAK-A applies the vector-sum Column- Cholesky algorithm to factor a general 
sparse matrix. Readers are referred to [12] for a comprehensive description of various storage 
schemes which result in efficient implementations of the algorithm. 

For n~ 5 and j — 3, the in-place Column- Cholesky factorization scheme thus transforms 
Oi t 3 to 4,3 for 3 < i < 5 as depicted in Figure 18. Note that the elements actually involved 
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Figure 18: Computing the third column of the Cholesky factor L . 
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in computing the third column of L, denoted as i* 3 , in the above example are shown in 
Figure 19. They are the coefficients of the third column of A and those of the computed L 
with their row indices greater than or equal to 3. Liu [15] makes the observation that if A is 
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Figure 19: The Uj's accessed and the a^j’s modified in computing L* 3 

read into memory one column at a time and each column of L is written out to the auxiliary 
storage as soon as it is computed, the in-core working space can be economized by keeping 
only those s which are needed for the current stage of elimination. Suppose the computed 
lij's are saved in a linear array sequentially, we use the above example to demonstrate the 
necessary data reorganization when the size of this working array is LNZSZE = 9. As 
shown in Figure 20, the lij elements are relocated (by overwriting elements which are not 
accessed any more) to make room for the newly computed Uj’s. For this example, data 
reorganization is necessary only before computing the third column and the fourth column 
of L . Clearly, the larger the size of the working array the fewer number of times the data 
reorganization needs to be done. 
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Figure 20: The organization of lijs in the working array. 
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In reference [15], Liu applies the idea above to large sparse matrices in his development of 
an adaptive general sparse out-of-core Cholesky factorization scheme. One of the advantages 
the algorithm features is that the frequency of data structure reorganization is adaptive to 
the available working space. Liu’s implementation of the out-of-core Cholesky scheme is 
compatible with SPARSPAK-A and is intended to be used in Step 3. We have incorporated 
this set of subroutines into an experimental processor in the Testbed and we shall report 
its performance on a set of CSM Testbed demonstration problems in Section 5. 

4.2,2 The Design of the Processor SPK 

Currently the entire interface together with the driver and a subset of SPARSPAK-A mod- 
ules are installed as a single processor SPK which can be invoked by the macroprocessor 
command [XQT SPK during the execution of the Testbed. The choice provided by this 
particular subset of SPARSPAK-A modules is the “Minimum Degree ordering [17]”. Fol- 
lowing the guideline contained in Section 6.2.1 of reference [22] for coding new processors, 
the main program of the processor SPK is implemented as a subroutine (named “5PA”) 
called by the Testbed executive module “ NICESPAR ”. Referring to the control diagram 
given in Figure 21, observe that the subroutine SPK calls another module “ SPKA ” which 
serves as the driver of SPARSPAK-A modules. In short, the subroutine SPKA allocates 
memory, sets up the problem by calling CSM-Interface modules, and solves the problem by 
calling SPARSPAK-A computational modules. The role the CSM-interface modules play 
is to retrieve the assembled linear system to be solved from the global database and input 
the problem in an appropriate form to SPARSPAK-A. The design of the processor at this 
level is thus generic and may be used with other sparse matrix packages. 

The CSM-interface consists of twenty-two modules. For easy reference, we list the 
subroutine or function name of each module and its formal parameters (if there is any) in 
Table 2 together with those of the two driver subroutines SPK and SPKA . All of these 
modules are written in the FORTRAN 77 language and a complete listing of programs 
is provided in Appendix C of this report. We shall discuss some implementation issues 
in section Section 4.2.3 and describe how these modules interface with the Testbed global 
database and SPARSPAK-A in Sections 4.2.4 and 4.2.5. The usage of the interface is 
described in Section 4.3. 


32 




Figure 21: The control diagram of the new processor SPK. 
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Table 2: The SPK driver and interface modules. 


DRIVERS 


SUBROUTINE SPK 
SUBROUTINE SPKA (A, MXSTOR) 
DOUBLE PRECISION A(l) 
INTEGER MXSTOR 


CSM-INTERFACE INITIALIZATION MODULES 

SUBROUTINE SPKCSM 

REAL FUNCTION CTIME ( IDUMMY) 

INTEGER IDUMMY 


CSM-INTERFACE PROBLEM INPUT MODULES 

SUBROUTINE GET JDF ( IBUF ) 

INTEGERS IBUF(l) 

SUBROUTINE GETDOF ( DOF, IBUF ) 

INTEGER*4 DOF(l), IBUF(l) 

SUBROUTINE GTZERO ( DOF, FBUF, MASK ) 

DOUBLE PRECISION FBUF(l) 

INTEGERS MASK(l), DOF(l) 

SUBROUTINE GTCOND ( DOF, IBUF, KC, MASK, CSIZE ) 
INTEGER*4 DOF(l), IBUFf 1 ), K C( 1 ), MASK(l), CSIZE 
SUBROUTINE GTMOTI ( FBUF, MASK, FCON, CSIZE ) 
INTEGER*4 MASK(l), CSIZE 

DOUBLE PRECISION FBUF(l), FCON(l) 

SUBROUTINE GETIJ ( DOF, IBUF, ICLQ, MASK, S ) 
INTBGER*4 DOF(l), IBUFfl), ICLQ(l), MASKfl), S(l) 
SUBROUTINE GTFORC ( FBUF, MASK, S ) 

INTEGERS MASK(l) 

DOUBLE PRECISION FBUF(l), S(I) 

SUBROUTINE GTNUM6 ( DOF, FBUF, MASK, FCON, S ) 
INTEGER*4 DOF(l), MASK(l) 

DOUBLE PRECISION FBUFfl), FCON(l), Sfl) 

CSM-INTERFACE UTILITY MODULES 

INTEGER FUNCTION SPACE ( IDUMMY ) 

INTEGERS IDUMMY 
SUBROUTINE LIBOPN 
SUBROUTINE QKINFO ( DSNAME) 

CHARACTER* 5 1 DSNAME 

SUBROUTINE GTRECI ( RECNUM, IBUF, LEN ) 
INTEGER*4 RECNUM, IBUF(l), LEN 
SUBROUTINE GTRECF ( RECNUM, FBUF, LEN ) 
INTEGER*4 RECNUM, LEN 

DOUBLE PRECISION FBUF(l) 

CSM-INTERFACE ERROR HANDLING MODULES 

SUBROUTINE EMSG 
SUBROUTINE EMSGO 
SUBROUTINE EMSG1 

SUBROUTINE EMSG2 

SUBROUTINE DEMSGO 


CSM-INTERFACE STATISTICS MODULES 


SUBROUTINE GETSOL ( FBUF, SOL, RATIO ) 
DOUBLE PRECISION FBUF(l), SOL(l), RATIO 
SUBROUTINE STATCS 
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4.2.3 Implementation Issues 

The two implementation issues we shall discuss in this section are “memory allocation” and 
“module/module communication”. 

Memory allocation Firstly, we note that the maximum working array storage available to 
the processor SPK is determined by the blank common dimension identically declared 
in the Testbed executive NICESPAR and the subroutine SPK , namely 

COMMON A(KSZZZ) 

Consequently, if the number of words provided by the blank common is insufficient 
for the processor SPK to solve a particular problem in-core, the dimension of the blank 
common must be increased, and the Testbed and the subroutine SPK must both be 
recompiled. 

We supply blank common of dimension KSZZZ (words) to the subroutine SPKA 
as a floating-point array of dimension MXSTOR . To accomplish this, we have the 
subroutine SPK execute the following statement: 

CALL SPKA ( A, MXSTOR ) 

where the value of MXSTOR is either KSZZZ or KSZZZ/2 depending on whether A 
is declared as a single-precision or double-precision array in the subroutine SPKA . 

All integer and floating-point arrays required by the CSM-Interface modules and 
SPARSPAK-A are then allocated by the subroutine SPKA from the one dimensional 
floating-point array A(MXSTOR). Note that in order to interact with SPARSPAK- 
A, the user is required to pass a working array 5 to the package and the location 
of 5 is the only parameter appearing in all of the SPARSPAK-A interface modules. 
In our case, the array S must be allocated from the working array A(MXSTOR). 
We have thus divided A(MXSTOR) into two segments. The top segment accom- 
modates arrays to be passed to the CSM-interface modules and the entire bottom 
segment is passed to SPARSPAK-A. If we let the variable MX USED denote the size 
of the top segment, the parameter to be passed to SPARSPAK-A is A (SPK), where 
SPK = MXUSED+ 1. 

A labeled common block CSMMAP is designated to keep the locations (origins in 
A) of the various arrays. The variables in COMMON /CSMMAP/ and the relative 
locations they represent are depicted in Figure 22. The type and size of the working 
arrays are tabulated in Table 3. Note that the buffer space for reading integer and 
floating-point records has been overlapped. 
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Table 3: The type and size of the SPK working arrays. 


Type 

Formal 

parameter 

Actual 

parameter 

Size 

Comments 

INTEGERS 

DOF 

A(DOF) 

NUMJNT + 1 

NUMJNT= total# of joints 


MASK 

A (MASK) 

NEQNS 

NEQNS = total# of equations 


KC 

A(KC) 

MAXDOF+1 

MAXDOF = 6 


ICLQ 

A (ICLQ) 

MAX DOF 



IBUF 

A(BUF) 

BUFMAX 

maximum buffer length 

DOUBLE PRECISION 

FBUF 

A(BUF) 

BUFMAX 

maximum buffer length 


FCON 

A(FCON) 

CSIZE 

total # of nonzero constraints 


SPK 

A (SPK) 

MAX ST OR— SPK -f 1 

the bottom segment of A 


Module/module communication The following labeled common blocks have been used 
to organize the communication between the SPK modules and the CSM Testbed mod- 
ules, between the SPK modules and the SPAESPAK-A modules, and among the mod- 
ules within the interface. 

1. COMMON /IANDO/ IIN, IOUTX. The two integer variables contain user input 
and output unit numbers assigned by the Testbed subroutine INTRO when the 
new processor begins execution. 

The /IANDO/ common appears in the SPK initialization subroutine SPKCSM 
and the SPARSPAK-A initialization subroutine SPRSPK. 

2. COMMON /SPAUSR/ MSGLVA, IERRA, MAXSA, NVARS. The /SPAUSR/ 
common allows user and/or processor SPK to communicate with SPARSPAK-A 
or vice versa. The meaning of the four integer variables are explained in Section 
4.3.2 and Section 4.3.3. 

The /SPAUSR/ common appears in the SPK subroutine SPKA which serves as 
the driver of SPARSPAK-A. 

3. The following common blocks are for communication among the SPK modules. 

COMMON /CSMSYS/ (6 variables) 

COMMON /CSMSPK/ (6 variables) 

COMMON /CSMUSR/ (11 variables) 

COMMON /CSMMAP/ (7 variables) 

COMMON /CSM CON/ (4 variables) 

COMMON /CSMDTA/ (8 variables) 

COMMON /PRBLEM/ (3 variables) 

The collection of related variables into a labeled common block avoids passing 
long parameter lists in the use of the subroutines and yet makes the coupling 
between modules easy to identify. Comments on the variables contained in these 
labeled commons are made at appropriate places throughout Sections 4.2.4, 4.2.5 
and 4.3. 
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Figure 22: Storage allocation of the SPK working arrays. 
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labeled commons are made at appropriate places throughout Sections 4.2.4, 4.2.5 
and 4.3. 

4.2.4 Interfacing with the Global Database 

There are eight modules in the interface which retrieve data from the global database and 
process them. The names of these subroutines are “ GETJDF ”, “ GETIJ ”, “ GTZERO ”, 

“ GTCOND ”, “ GTFORC ”, “GTMOTP\ “ GTNUM5 ” and “ GETSOL ”. We shall use “Gxxxxx” 
to represent an arbitrary one of them. All of these modules retrieve datasets from the 
Testbed via two utility modules which are either “ QKINFO and GTRECI ” (for retrieving 
integer records) or “ QKINFO and GTRECF ” (for retrieving records containing floating- 
point numbers). Figure 23 depicts the coupling of the interface modules with the GAL- 
processors. Readers are referred to reference [26] for a complete description of the calling 
sequence and the operation of each GAL-processor employed. 



L . I L I L I I I I I I I 


( | Global 

database 

Dataset name 


Figure 23: The coupling of CSM-interface modules and GAL-processors. 
For each designated dataset, the labeled common /CSMSPK/ is used to 
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1. provide the input arguments LDI and TRACE to the GAL-processors. (The meaning 
of LDI and TRACE is given in Table 4.) 

2. store the dataset attributes the interface module QKINFO acquires from the GAL- 
processors LMFIND, GMEGKA and GMGECY. 

3. communicate the dataset attributes to the interface modules Gxxxx, and the GAL- 
processors GMCORN and GMGETN via the interface module GTRECI or GTRECF. 

The /CSMSPK/ common thus appears in QKINFO , GTRECI, GTRECI and each Gxxxx 
module. The variables contained in /CSMSPK/ and their meaning are given in Table 4. 


Table 4: The variables in COMMON /CSMSPK/. 


| COMMON /CSMSPK/ 

variable 

meaning 

IDSN 

LDI 

NLEN 

NREC 

RTYPE 

TRACE 

Dataset sequence number. 

Logical Device Index of library device. 

The record length. 

The number of records in the dataset. 

The data type. 

A positive integer used as identifying label 
in error traceback prints. 


Since the actual data contained in each dataset is unique, each subroutine Gxxxxx must 
be specifically coded to interpret the data retrieved. The datasets to be accessed by the 
eight interface modules are listed in Table 5. For each dataset, given in Table 5 are also 
the name of its source processor and the name of the dedicated interface module. The 
last column of Table 5 indicates the appropriate utility module which should be called to 
retrieve the type of data provided by the specified dataset. 

The data retrieved from each dataset and how they are handled by the interface routines 
are described below. Readers are referred to reference [22] for a description of the format 
of each dataset. 

JDF1.BTAB.1.8 provides the total number of joints and the maximum number of active 
(unconstrained) degrees of freedom a joint may have in the model. 

The subroutine GETJDF retrieves the data and stores them in the variables NUMJNT 
and MAXDOF in the labeled common 

/PRBLEM/ MAXDOF, NEQNS, NUMJNT 

In an attempt to be flexible in handling the more general case in the future, the sub- 
routine GETDOF stores the active degrees of freedom for each individual joint in an 
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Table 5: Datasets accessed by Gxxxxx and GTRECx. 


Source Processor 

Dataset 

Gxxxxx 

GTRECx 

TAB 

JDFl.BTAB. 1.8 

GETJDF 

GTRECI 

K 

K.SPAR.jdfB 

GTZERO 

GTRECF 

TAB 

CON. .neon 

GTCOND 

GTRECF 

T0P0 

KMAP..nsubs.ksize 

GETIJ 

GTRECI 

AUS 

A PPL. FOR C. iset. 1 

GTFORC 

GTRECF 

AUS 

A PPL. MO T I. iset. 1 

GTMOTI 

GTRECF 

K 

K. SPAR, jdf 2 

GTNUM5 

GTRECF 

SSOL 

S TA T. DISP. iset. neon 

GETSOL 

GTRECF 


accumulated form in an integer array DOF so that the number of degrees for joint 
number I can be computed from DOF(I+l)—DOF(I), where D0F(1)—1 , and that 
DOF(NUMJNT+l)—DOF(l) gives the total number of equations of the assembled 
system. The latter value is also stored in the variable NEQNS in the /PRBLEM/ com- 
mon. Since the current version of the CSM Testbed assumes uniform degrees of free- 
dom per joint in storing the system stiffness matrix, DOF(I+l) — DOF(I)=MAXDOF 
for 1 <1 <NUMJNT. 

K.SPAR.jdf2 provides the assembled global stiffness matrix stored as an array of JDFxJDF 
submatrices, where JDF is the maximum degrees of freedom in the model and its value 
is available from the the variable MAXDOF in the /PRBLEM/ common block. Note 
that the integer jdf 2 in the name of this dataset is the square of the value of JDF . 

The subroutine GTZERO retrieves the system stiffness matrix and detects dummy 
rows by examining its diagonal elements. For each zero diagonal coefficient detected, 
a zero is entered into the integer array MASK at the location MASK(I), where I 
is the equation number of the dummy row. The convention we have adopted is 
that MASK(J)— —1 if the J th equation is neither constrained nor a dummy row, 
MASK(J)— 0 if it corresponds to a dummy row or a zero constraint, MASK(J)~ 1 if 
it corresponds to a nonzero constraint. 

CON. .neon provides constraint information for each joint degree of freedom. The informa- 
tion available indicates for each joint which component is free, which component is 
constrained to be zero and which component has a non-zero constraint. Such informa- 
tion is encoded so that one integer is stored for each joint in the model. The current 
encoding mechanism assumes that the maximum number of degrees of freedom a joint 
may have is “six”. The constraints corresponding to the six degrees of freedom are 
encoded into the right most six bits of a seven-bit integer. The subroutine DECODE 
accepts an integer as input and returns the status of each of the MAXDOF degrees 
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of freedom in the leading MAXDOF locations of a working array of length seven. 

The subroutine GTCOND retrieves the encoded data from CON.. neon, calls DECODE 
to obtain the constraint status for each joint in the model, and sets the corresponding 
entries in the integer array MASK to be “0” or “1” as explained above. An inte- 
ger output parameter CSIZE records the total number of nonzero constraints whose 
numerical values are expected to be available in the dataset APPL.MOTI.iset.l . 

Therefore, after both subroutines GTZERO and GTCOND are executed, all con- 
straint information is available for other SPK modules in the integer array MASK. 
Note that we have treated the dummy rows as if they correspond to zero constraints. 

KMAP..nsubs.ksize provides the block nonzero structure of the system stiffness matrix. 
Note that the value of nsubs in the name of the dataset represents the total number 
of submatrices in the system stiffness matrix for the model, and that the integer ksize 
is the maximum number of joints active at any time during the assembly of the system 
matrix. 

The subroutine GETIJ accesses KMAP..nsubs.ksize and the integer array MASK to 
obtain the matrix structure for input to SPARSPAK-A. We explain how the con- 
straints are handled in Section 4.2.5. 

APPL.FORC.iset.l provides applied forces and moments on each joint in each active direc- 
tion. The integer iset in the dataset name identifies a unique load case. 

The subroutine GTFORC retrieves the data but inputs a retrieved numerical value 
as a component of the right-hand side vector to SPARSPAK-A only if it does not 
correspond to a variable constrained to be zero (i.e., MASK(I)^ 0 if / is the equation 
number). 

Since the right-hand side is initialized to be identically zero in SPARSPAK-A, and the 
modifications to the right-hand side caused by nonzero constraints are to be “added” 
to the appropriate components by subroutine GTNUM5 , the input of right-hand side 
to SPARSPAK-A is not completed before the subroutine GTNUM5 is executed. 

APPL.MOTI.iset.l provides applied motions on each joint in each active direction. As 
mentioned earlier, the integer neon in the name of this dataset identifies a particular 
constraint case, and numerical values for the nonzero constraints detected by the 
subroutine GTCOND are expected from this data set. 

The subroutine GTMOTI retrieves the available applied motions and stores them 
in a floating-point array FCON(I), where 1 < I < CSIZE, and CSIZE is the total 
number of nonzero constraints determined in the subroutine GTCOND. Therefore, 
when CSIZE = 0, the subroutine GTMOTI will return without attempting to access 
the dataset. However, when CSIZE> 0, it is a fatal error if the dataset is missing or 
less than CSIZE values are available. 
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STAT.DISP.isetncon provides the computed static displacements for each joint in each 
active direction. Unique solution is obtained by specifying the load set and constraint 
case in the name of the dataset. 

The subroutine GETSOL retrieves the Testbed solution from this dataset and verifies 
the correctness of the SPARSPAK-A solution by computing its relative error with 
respect to the Testbed solution. More details in this aspect are provided in Section 5 
on numerical experiments. 

4.2.5 Interfacing with SPARSPAK-A 

The processor SPK may interact with SPARSPAK-A via the interface modules given in 
Table 6, which correspond to our choice of the minimum degree ordering ( subroutine OR - 
DRB5) for the new processor. 


Table 6: SPARSPAK-A interface modules - a subset. 


Initialization of SPARSPAK-A 

SPRSPK 

Structure input 

IJBEGIN 
INIJ ( I, J, S) 

INROW ( I, NIR, IR, S) 
INIJIJ ( NIJ, II, JJ, S) 
INCLQ ( NCLQ, CLQ, S) 
IJEND( S ) 

Ordering 

0RDRB5 ( S ) 

Matrix input 

INAIJ5 ( I, J, VALUE, S) 

IN ROW 5 ( I, NIR, IR, VALUES, S) 
INMAT 5 ( NIJ, II, JJ, VALUES, S) 

Right-hand side input 

INBI(I, VALUE, S) 

INBIBI ( NI, II, VALUES, S) 
INRHS (RHS, S ) 

Factorization and/or Solution 

SOLVE5 ( S ) 

Relative error estimation 

EREST5 ( RELERR, S ) 

Print statistics 

STATSA 

Save and Restart the computation 

SAVE A (K, S) 
RSTRTA ( K, S) 
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The coupling of the SPK modules and SPARSPAK-A is depicted in Figure 24. The mod- 
ules which interact with SPARSPAK-A are “SPAA”, “GETIJ”, “ GTFORC ” and “GT- 
NUM5 ”. The module SPKA serves as the driver program of SPARSPAK-A. The module 
GETIJ inputs the nonzero structure of the system stiffness matrix to SPARSPAK-A. The 
modules GTFORC and GTNUM5 are involved in inputting nonzero coefficients and the 
right-hand side to SPARSPAK-A. The particular SPARSPAK-A subroutines to be called 
by each of these interface modules are explicitly given inside the dashed boxes. 
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Figure 24: The coupling of the processor SPK and SPARSPAK-A. 
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Since SPAESPAK-A modules do not handle constraints, the retrieved system stiffness 
matrix and the right-hand side must be adjusted before they can be input to SPARSPAK-A. 
The necessary modifications to the structure and the numerical values are detailed below. 

Input the structure of the system stiffness matrix to SPARSPAK-A - In this sec- 
tion we describe how the subroutine GETIJ inputs the the structure of the system 
stiffness matrix to SPARSPAK-A. The dataset KMAP..nsubs.ksize contains the sys- 
tem topology map. From this map we can retrieve the following information for each 
joint. 

JNT - The number of the current joint. 

CONRNG - The number of submatrices including the diagonal in the upper triangle 
for the current joint. 

CONECT(CONRNG-l) - A list of joints connected to the current joint. 

Let us consider the following finite- element model which is given as an example in 
reference [23]. 


5 6 


12 3 4 

Figure 25: A model. 


Table 7: A model. 


Element 
# type 

Connected 

Nodes 

1 BEAM 

1,2 

2 

2,3 

3 

3,4 

4 

2,5 

5 

3,6 

6 PLATE 

1,2,5 

7 

2, 3, 6, 5 

8 

3,4,6 


For this model, the information expected to be available in KMAP..nsubs.ksize is 
listed in Table 8. 
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Table 8: From dataset KMAP..nsubs.ksize. 


JNT 

CONRNG 

CONECT(CONRNG-l) 

1 

3 

2,5 

2 

4 

3, 5, 6 

3 

4 

4,5, 6 

4 

2 

6 

5 

2 

6 

6 

1 



Given in Figure 26 is the upper triangular block structure of the system matrix (in- 
cluding the diagonal blocks) described by Table 8. 
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Figure 26: Upper triangular block structure of the system matrix for the model 


problem. 
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If each joint has three degrees of freedom in the model, each block is a 3 X 3 submatrix 
and the system stiffness matrix K has the nonzero entries as given in Figure 27. 
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Figure 27: Nonzero entries in the upper triangle of K (including diagonal submatrices.) 


If every degree of freedom is active (unconstrained) on each joint, then the structure 
input to SPARSPAK-A is as specified in Figure 27. It should be pointed out that 
because SPARSPAK-A anticipates only “symmetric” nonzero structure, the struc- 
ture input routine always records a logical nonzero in both (i,j) and (ji, i) positions 
regardless of which index pair is actually being entered. Furthermore, the package 
automatically removes duplications so that it does not matter if both (z, j) and (jf, i) 
pairs are entered. 

In order to demonstrate how we handle the constrained degrees of freedom, let us 
assume that the second degree of freedom on joint number 5 is constrained. In this 
case, the corresponding columns and rows of data in K except for the diagonal elements 
will be treated as zero entries. The nonzero positions SPARSPAK-A is informed of 
consist of the remaining nonzeros as given in Figure 28. 
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Figure 28: Remaining nonzero entries in the upper triangle of K. 


As seen from Figure 28, the equations corresponding to the constrained degree of freedom 
is the fourteenth equation. We have thus ignored the nonzero entries in locations (i, 14) 
and (14, i) for all V s except for the diagonal entries. Accordingly, the numerical coefficients 
corresponding to these ignored locations must not be input to SPARSPAK-A and the right- 
hand side must he appropriately adjusted to reflect the change of the system matrix. We 
next explain the internal working of our numerical input module. 

Input the numerical values to SPARSPAK-A - The subroutine which inputs the nu- 
merical values to SPARSPAK-A and modifies the right-hand side according to each 
constrained degree of freedom is GTNUMi , where i = 1,3, and 5 distinguishes the 
SPARSPAK input modules INAIJi called for each ordering. 

To see how the right-hand side should be modified, we refer to Figure 29 for the same 
example, where we label each ignored coefficient a t -j explicitly, and indicate that the 
coefficient for the diagonal entry ai 4 i 4 is set to 1. 

Let the nonzero constraint corresponding to the second degree of freedom on joint 
number 5 be ci 4 . Our change to the system matrix and right-hand side should reflect 
the following. 

1. The fourteenth equation is replaced by 


Zl4 = c x4 . 
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Figure 29: Nonzero entries in the upper triangle of K. 
2. Modify the right-hand side to be 



Thus, the right-hand side elements i — 1,2, 3, 4, 5, 6, 7, 8, 9, 13 are modified to 
be 

^ X Cj4 
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and the right-hand side element bj , j = 15, 16, 17, 18 must be modified to be 

* bj — Ui4j X Ci 4 . 

To summarize, for each coefficient retrieved from the dataset K.SPAR..* , sub- 
routine GTNUMi checks whether the corresponding degree of freedom is constrained. 
If that is the case, the value of 1.0 will be input to SPARSPAK-A as a^i and the 
constraint value is input to SPARSPAK-A as b{. 

For each off-diagonal element GTNUMi checks the following four possible cases. 

1. If both and Xj are constrained, no action needs to be taken. 

2. If X{ is active and Xj is constrained to be cj then 

b{ < b{ — d{ j X Cj . 

3. If X{ is constrained to be C{ and Xj is active then modify 

bj 4 bj ' cijj X C{ . 

4. If neither X{ nor Xj is constrained, input the retrieved value to SPARSPAK-A 

and specify the location to be (j, i). (SPARSPAK-A requires the numerical value 
to be input for the lower triangular part only.) 
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4.3 The Usage of the Interface 
4.3.1 The Execution Path 

The usage of the interface in solving a Testbed problem is reflected by the execution path of 
the subroutine SPKA as sketched in Table 9. The execution sequence is enforced by checking 
and updating the value of the variable STAGE in the common block /CSMCON/. The 
values of STAGE for the successful completion of each corresponding step are listed in the 
last column of Table 9. 


Table 9: The execution path of the subroutine SPKA 


1 Execution path 1 

SPARSPAK-A 

Interface 

Dataset 

/CSUCONf 



subroutine 

subroutine 

dependency 

STAGE 

Step 1.1 

Start 

SPRSPK 

SPKCSM 


0 

1.2 



LIB OP N 


10 

1.3 



CTIME 



Step 2.1 

Problem input 


GETJDF 

JDFl.BTAB.1.8 

20 

2.2 



GETDOF 


30 

2.3 



GTZERO 

K.SPAR.jdf2 

40 

2.4 



GTCOND 

CON. .neon 

50 

2.5 



GTMOTI 

APPL.MOTI.ittU 

60 

Step 3 

Structure input 

IJBEGIN 

GETIJ 

KMAP..nsubs.ksize 

70 



INCLQ 






INIJ 






IJEND 




Step 4 

Order and 

ORDRB5 





allocate storage 





Step 5 

Input numerical 

INBI 

GTFORC 

APPL.FORC.iset.l 

80 


values for b 





Step 6 

Input numerical 

1NAIJ5 

GTNUM5 

K.SPAR.jdf2 

90 


values for A and b 

INBI 




Step 7 

Factor A and solve 

SOLVE5 





for solution x 





Step 9 

Relative error estimation 

E REST 5 




(optional) 






Step 10 

Compare x with 


GETSOL 

STAT.DISP.isei.ncon 


(optional) 

CSM Testbed solution 





Step 11 

Collect statistics 

STATS A 

STATCS 



(optional) 
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4.3.2 User Input to the Processor SPK 

In our current implementation of the processor SPK, the user-processor communication is 
accomplished using an external text file. The input requirement and format are reflected 
by the following code segment of the subroutine SPKA 

c 

SUBROUTINE SPRA ( A, MXSTOR ) 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

12 

C 


INDADA = 41 


SET MSGLVL AS DESIRED 


READ ( INDATA, 12 ) MSGLVL 


SET MSGLVA AS DESIRED 


READ ( INDATA, 12 ) MSGLVA 


SET MAXIMUM BUFFER LENGTH 


READ ( INDATA, 12 ) BUFMAX 
FORMAT ( 14 ) 


INPUT NAME OF LIBRARY AND DATASETS FOR GIVEN PROBLEM 


READ 

( 

INDATA, 

22 

) 

LIBNAM 

FORMAT ( A40 ) 




READ 

( INDATA, 

32 

) 

JDFSET 

READ 

( 

INDATA, 

32 

) 

KMAP 

READ 

( 

INDATA, 

32 

) 

KSPAR 

READ 

( 

INDATA, 

32 

) 

CON 

READ 

( 

INDATA, 

32 

) 

APPLF 

READ 

( INDATA, 

32 

) 

APPLM 

READ 

( 

INDATA, 

32 

) 

STATD 


32 FORMAT ( AS1 ) 


RETURN 

C 

END 

The following comments are in order. 

1. As shown in the above code segment, we have designated the logical unit number 
41 to be used for the input data file. This choice is made under the restriction that 
logical unit numbers 1 through 40 should not be used for files other than libraries to 
avoid possible conflicts with CLIP and GAL [22]. 
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2. The variable MSGLVA stands for “message level of SPARSPAK-A”. The user may 
govern the amount of output from SPARSPAK-A by setting MSGLVA to the values 
Table 10. 


Table 10: The valid input values of MSGLVA . 


MSGLVA 

amount of output 

0 

no information is provided. 

1 

only warnings and errors are printed. 

2 

warnings, errors and summary are printed. 

3 

warnings, errors, summary and some statistics are printed. 

4 

detailed information for debugging purposes. 


3. The variable MSGLVL allows user to control the amount of output from the interface 
modules. Given in Table 11 are the input values acceptable for MSGLVL . 


Table 11: The valid input values of MSGLVL . 


MSGLVL 

amount of output 

0,1 

no information is provided. 

2 

warnings, errors and summary are printed. 

3 

detailed information for debugging purposes. 


4. The value of BUFMAX should be set to the maximum record length of any dataset 
the processor SPK ever needs to retrieve. 

5. The variables initialized by user input are collected into the two labeled common 
/SPAUSR/ and /CSMUSR/. 

6. An example - To solve the linear system of the test problem demol using SPARSPAK- 
A, edit a file named “fort.41” to contain the following data: 

2 

2 

2240 

/uar . MC 6 8020/nlal / echu/ns/DEMO/demo 1.101 

JDF1.BTAB.1.8 

KMAP. .9.3 

K. SPAR. 36 

C0N..1 

APPL.F0RC.1 .1 
APPL.MOTI. 1 . 1 
STAT.DISP.l . 1 
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Note that the path name of the library file “demol.101” is installation dependent. The 
dataset names listed above can be identified from the table of contents of the library 
demol .101 given in Figure 3. Note that the datasets APPL.FORC.iset.l and APPL.MOTLisetl 
may not both exist, and it is indeed the case for the problem demol - one cannot find the 
name APPL.MOTI.1.1 listed in the table of contents of its data library. However, as noted 
above, we have required the user to input both names in order to maintain a uniform format 
for user input. In this case, the variable APPLM is simply a dummy variable, because the 
subroutine GTMOTI will not attempt to access this dataset as explained in Section 4.2.4. 

4.3.3 Output from the Processor SPK 

1. Output from SPARSPAK-A: Readers are referred to Section 7 of the SPARSPAK-A 
User’s Guide [2] for a complete description of the statistics and error messages output. 

2. Output from the interface modules: 

(a) Statistics gathering (STATCS) - The information contained in Table 12 may be 
printed by the following statement. 

CALL STATCS 


Table 12: Information printed by the subroutine STATCS. 


MSGLVL 

Information 

Variable 

Common block 

1 

Total CSM-time required 
Maximum CSM-storage required 

CSMTIM 

CSMSTR 

/CSMDTA/ 

2,3 

Size of storage array 

MAXCSM 

/CSMUSR/ | 

2,3 

Number of joints 

Max degree of freedom per joint 

Number of equations 

NUMJNT 

MAXDOF 

NEQNS 

/PRBLEM/ 

3 

Addresses of arrays 

DOF 

BUF 

MASK 

KC 

ICLQ 

FCON 

SPK 

/CSMMAP/ 


(b) Error messages ( IERR ) - When fatal error is detected, so that the computation 
cannot proceed, a positive code is assigned to the variable IERR in the common 
block /CSMUSR/ . The names of the modules in which the error occurs, the 
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Table 13: Error messages of the processor SPK. 



IERR 

Error message 

SPACE 

1001 

Insufficient storage. The last stage completed 
and the required storage are printed 

LIBOPN 

MjBM 

Cannot open dataset library 


Eh 

The maximum logical device index = 30. The LDI 
returned exceeds this value. 

GETJDF 

■SSI 

Incorrect execution sequence. 


■El 

Dataset does not have all expected items. 

GETDOF 


Incorrect execution sequence. 

GTZERO 

I 

Incorrect execution sequence. 

GTCOND 

I 

Incorrect execution sequence. 

GETIJ 

1023 

Incorrect execution sequence. 

GTFORC 

1024 

Incorrect execution sequence. 

GTMOTI 


Incorrect execution sequence. 



Unexpected nonzero constraint value. 



Zero entry for a nonzero constraint occurs. 

GTNUM5 

1028 

Incorrect execution sequence. 

QKINFO 


LMFIND : cannot find dataset. 



GMGEKA: record does not exist. 


1 

GMGECY : record group key undefined. 



GMGECY : segmented record group noted. 


2009 

Insufficient buffer space. The required value 
for the input variable B UFMAX is printed 

GTRECI 

ESSM 

record type mismatch • • • 


EH 

GMGETN: error detected by LMERCD ••• 

wsssmBm 


record type mismatch ■ • • 

■Hill 


GMGETN: error detected by LMERCD ■ ■ ■ 


numerical error codes, and the corresponding error messages as given in Table 13 
may be printed by setting the variable MSGLVL to be “2” or a higher number. 
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4.3.4 An Example - Solving the Testbed problem demol 
Input data: 

2 

2 

2240 

/usr. MC68020 /nlal / echu/ns / DEMO/ demol . 101 

JDF1 . BTAB . 1 .8 

KMAP..9.3 

X. SPAR. 36 

COH. .1 

APPL.FORC.1.1 

APPL.MOTI.1.1 

STAT.DISP.1.1 

The following output is produced by the macroprocessor command [xqt SPK 

** BEGIN SPK ** DATA SPACE= 600000 WORDS 

1 

********** UNIVERSITY OF WATERLOO 
********** SPARSE MATRIX PACKAGE 
********** (SPARSPAK) 

********** RELEASE 3 

********** (c) JANUARY 1984 

********** ANSI FORTRAN 
********** DOUBLE PRECISION 
********** LAST UPDATE JANUARY 1984 


OUTPUT UNIT FOR ERROR MESSAGES 6 

OUTPUT UNIT FOR STATISTICS 6 

LIBOPI- OPEN /usr .MC68020/nlal/echu/ns/DEM0/demol .101 
<DM> OPEN, Ldi: 2, File: /u s r. MC6802 0 /nlal /echu/ns/ DEMO /demol .101 , 

Attr: rold. Block I/O 

DATASETS TO BE ACCESSED: 

JDF1. BTAB. 1.8 
KMAP. .9.3 
K. SPAR. 36 
CON. .1 

APPL.FORC.1.1 

APPL.MOTI.1.1 

STAT.DISP.1.1 

GETJDF - GET NUMBER OF JOINTS AND . . . 

GETDOF - GET DEGREES OF FREEDOM . . . 

GTZERO - DETECT DUMMY ROWS . . . 

GTCOND - GET CONSTRAINTED VARIABLES... 

GTMOTI - GET NONZERO CONSTRAINTS... 
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GETIJ - INPUT NONZERO STRUCTURES... 

IJBEGN - BEGIN STRUCTURE INPUT . . . 

INIJ - INPUT OF ADJACENCY PAIRS . . . 

IJEND - END OF STRUCTURE INPUT . . . 

ORDRBS - MINIMUM DEGREE ORDERING . . . 

GTFORC - INPUT RIGHT HAND SIDE... 

INBI - INPUT OF RIGHT HAND SIDE . . . 

GTNUM5 - GET NONZERO NUMERIC... 

INAIJ5 - INPUT OF MATRIX COMPONENTS . . . 

SOLVES - GENERAL SPARSE SOLVE ... 

ERESTS - ERROR ESTIMATOR . . . 

GETSOL - COMPARE WITH TESTBED SOLN . . . 

MAX. REL ERR COMPARED TO STAT. DISP. 1 . 1 
IS 0.4824782«-0T IN COMPONENT 26 

CSM SOL a 0.285208672285080+00 WE HAVE 0.285208686045780+00 
STATCS - SYSTEM-CSM STATISTICS . . . 


SIZE OF STORAGE ARRAY (MAXCSM) 300000 

NUMBER OF JOINTS 5 

MAX DEGREES OF FREEDOM PER JOINT 6 

NUMBER OF EQUATIONS 30 

TOTAL CSM-TIME REQUIRED 3.740 

MAXIMUM CSM-STORAGE REQUIRED 2271. 

STATSA - SYSTEM-A STATISTICS ... 

SIZE OF STORAGE ARRAY (MAXSA) 297729 

NUMBER OF EQUATIONS 30 

NUMBER OF OFF-DIAGONAL NONZEROS 336 

TIME FOR ORDERING 0.020 

STORAGE FOR ORDERING 442. 

TIME FOR ALLOCATION 0.000 

STORAGE FOR ALLOCATION 308. 

STORAGE FOR SOLUTION 367. 

TIME FOR FACTORIZATION 0.040 

TIME FOR SOLUTION 0.020 

OPERATIONS IN FACTORIZATION 956. 

OPERATIONS IN SOLUTION 396. 

TIME FOR ESTIMATING RELATIVE ERROR 0.040 

OPERATIONS IN ESTIMATING REL ERROR 1330. 

STORAGE FOR ESTIMATING REL ERROR 397. 

ESTIMATE OF RELATIVE ERROR 2.O880-O8 


57 


TOTAL TIME REQUIRED 
MAXIMUM STORAGE REQUIRED 
EXIT SPK CPUTIME= 4.2 I/0(DIR,BUF)= 


0 

442 

0 0 


120 
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5 Numerical Experiments 

In this section, we report experimental results of several matrix factorization processors we 
have installed in the CSM Testbed. 

5.1 The Specifications of the Test Problems 

For all processors, the tests are performed on the NICE/SPAR demonstration problems 
listed in Table 14. The finite element model of CSM focus problem 1 has been refined to 
generate larger problems focusl, focus2, focus3 and focus4. The five different meshes we 
have used are given in Table 15. 


Table 14: NICE/ SPAR demonstration problems. 


NICE/ SPAR demonstration problems 

p648 

CSM focus problem 1 — Buckling of a blade-stiffened 
panel with a discontinuous stiffener 

focusl 

p648 with finer mesh I 

focus2 

p648 with finer mesh II 

focus3 

p648 with finer mesh III 

focus4 

p648 with finer mesh IV 

demol 

Beam problem 

demo 2 

Vibration of a circular membrane 

demo3 

Circular plate problems 

demo4 

Rectangular plate problems 

demo6 

Cylindrical shell problems 

demo 7 

Buckling of a cylindrical shell due to torsional loading 

demo9 

Beam problems 

demo 10 

Saturn 5 Launcher Umbilical Tower (LUT) 

demo 12 

Hyperbolic paraboloid static solution 

demol 3 

Composite toroidal shell 


Each problem is completely specified by the datasets in Table 16 except that the load 
set APPL.FORC.iset.1 and the applied displacement dataset APPL.MOTI.iset.l may not 
both exist. For example, there is no applied force vector for the panel focus problem and 
there is no applied displacements for the static analysis of the mast problem. The value of 
neon selects one of possibly more than one constraint cases and the value of iset specifies a 
particular load case of applied force and moments, which is also the load case of the applied 
motions if there exist nonzero constraints. Corresponding to each pair of (neon, iset) there 
is a unique solution which may be retrieved from the dataset STAT .DISP .iset. neon to verify 
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Table 15: User-specified meshes for CSM focus problem 1. 


User-specified meshes for CSM focus problem 1 


NRINGS 

NSPOKES 

NELX 



NELS 

p648 

2 

8 

3 

1 

1 

1 

mesh I 

4 

16 

3 

1 

1 

1 

mesh II 

2 

8 

6 

2 

2 

2 

mesh III 

2 

8 

12 

2 

2 

2 

mesh IV 

4 

16 

6 

2 

2 

2 


the correctness of an experimental processor. The full names of the datasets can be found in 
the table of contents of the data library which can be looked up during or after the execution 
of a particular analysis in the Testbed. As shown in the example given in Table 16, a “0” 
component in the dataset name can be represented by a null entry. A sample content list of 
the data library demol .101 is given in Section 2 of this report, which was produced by the 
CLAMP directive *T0C during the execution of problem demol. For each test problem, the 
path name of its data library and the names of the datasets in Table 16 consist of the user 
input to an experimental processor. Note that the use of * as a component of the dataset 
name implies a generic wild-card match, hence it should not be used unless the dataset with 
its name matching the remaining components is unique in the data library. 

Table 16: Data sets accessed by CSM-SPARSPAK interface modifies. 


The accessed CSM Testbed datasets 

Name 

An example 
(neon, iset) = (3,6) 

JDF1.BTAB.1.8 

JDF1.BTAB.* 

KMAP. . nsubs.ksize 

KMAP..* 

K.SPAR.jdf2 

K.SPAR.* 

CON. .neon 

CON.. 3 

A PPL. FOR C. iset. 1 

CAPPL.FORC.6.1 

A PPL. MO TI. iset. 1 

APPL.MOTI.6.1 

S TA T. D1SP. iset. neon 

STAT.DISP.6.3 


The system Ax ~ b presented to each experimental processor is the upper triangular 
part of the system stiffness matrix retrieved from the dataset K.SPAR.jdf2 subject to the 
changes necessitated by the way we handle constraints and dummy rows. The modified 
system has the following characteristics. 

1. The coefficient matrix and the right-hand side are modified so that each equation 
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corresponding to a constrained variable X{ can be replaced by 

x i = c i > 

where c; > 0 is the specified constraint. 

2. The identically zero rows are detected before problem input and the corresponding 
variables are treated as being constrained to zero. 

3. The dimension of the modified coefficient matrix is equal to the product of the number 
of joints and the degree of freedom per joint in the model. The number of equations of 
each demonstration problem is given in Table 17 under the column heading “neqns”. 

In Table 17, we summarize the characteristics of the linear systems retrieved for each 
demonstration problem. The entries in the column labeled nonzeros in K.SPAR ” are 
computed from nsubs X jdf 2, where we recall that nsubs is the total number of submatrices 
in the block upper triangular part (including the diagonal blocks) of the system stiffness 
matrix and that jdf 2 =JDFxJDF represents the number of elements in each submatrix. 
Therefore, the nonzero count here includes the coefficients in the lower triangular part of the 
diagonal blocks and the coefficients in the dummy rows as well as the rows corresponding 
to the constrained variables. The actual off-diagonal nonzero elements input to an experi- 
mental processor are listed in the last column under the heading of off-diag nonzeros in 
A”. 


Table 17: Characteristics of NICE/SPAR demonstration systems. 


| Characteristics of the linear systems Ax = 6 ,j 



d.o.f 

neqns 

# dummy 
rows 

# zero 
constraints 

# nonzero 
constraints 

# nonzeros 
in K.SPAR 

# off-diag 
nonzeros in A 

I p648 

108 

6 

648 

78 

98 

10 

17064 

9706 

iesssbi 

192 

6 

1152 

154 

116 

12 

31320 

18458 | 


276 

6 

1656 

228 

167 

17 

45792 

26824 

focus3 

480 

6 

2880 

396 

167 

17 

81216 

50560 

focus4 

388 

6 

2328 

332 

185 

19 

65088 

38656 

demol 

5 

6 

30 

0 

6 

0 

324 

168 


101 

3 

303 

101 

203 

0 

4077 

342 

demo 3 

101 

5 

505 

0 

80 

0 

11325 

7830 

demo 4 

54 

5 

270 

0 

55 

0 

5675 

3546 

demo 6 

121 

6 

726 

0 

97 

0 

19476 

14151 

demo 7 

132 

6 

792 

0 

36 

0 

22464 

18576 | 

demo 9 

11 

6 

66 

0 

6 

0 

756 

474 


372 

6 

2232 

0 

24 

0 

47376 

39072 

demo 12 

36 

6 


0 

18 

0 

5256 

3978 

demo 1 3 

337 

6 


0 

96 

0 

59364 

49743 | 
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5.2 The Numerical Properties of the Test Problems 
5.2.1 The Conditioning of the System Stiffness Matrix 

In Table 18, we list the estimated condition number of the system stiffness matrix for each 
test problem. The condition numbers are provided by SPARSPAK-A and their computation 
is described in reference [1]. The order of magnitude of the condition numbers indicates 
that the single-precision solution of these problems my not have significant digits in some 
components. By comparing the single-precision static displacement solutions obtained from 
the Testbed processors INV and SSOL for the same problem using different joint orderings, 
our numerical experiments confirm that the loss of all significant digits can indeed occur in 
small components of the solution. 

Table 18: Numerical properties of NICE/SPAR demonstration problems. 


Condition number of 
the system stiffness matrix 


SPARSPAK-A estimate of 
condition number 

p648 

2.2 x 10 7 


3.7 x 10 7 ' 



focus3 


focus4 


demol 


demo2 

2.2 X 10 7 

demo3 

1.7 X 10 7 

lEZZZZM 

1.8 x 10 7 1 


qq^hii 

demo7 

3.2 x 10 7 

demo9 


demol 0 


demo 12 

5.6 x 10 9 

demo 13 

1.4 x 10 7 f 


5.2.2 The Accuracy of the Computed Solutions 

The condition number estimates we presented in Table 18 indicate that in order to have sig- 
nificant digits in all components of the solution to be stored in the dataset STAT.DISP.iseLncon, 
the system stiffness matrix should be stored in double-precision and processors INV and SSOL 
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should perform the numerical computation in double-precision. The following information 
from reference [24] tells us how to ensure that the computations by each processor are 
performed with the desired precision. 

1. Processor K stores the system stiffness matrix in double precision if the user input 
parameter SPDP is reset to 2 as shown in a sample script given later in this paragraph. 

2. Processor INV calculates the triangular matrix using double precision if the input 
system stiffness matrix dataset is stored in double precision. However, the factors 
output by processor INV will be stored in the precision determined by resetting the 
user- controlled argument SPDP: 1 (default) = single precision, 2 = double precision. 

3. Processor SSOL computes the displacement solution vector in double-precision if the 
factored matrix is stored in double-precision. 

Therefore, each reset SPDP in the following script excerpt ensures that the output dataset 
is in double precision, which in turn ensures that the computation by the next processor is 
performed in double precision. 


[xqt K 

reset SPDP=2 
[xqt INV 

reset SPDP=2 
[xqt SSOL 


For each demonstration problem, the solution provided by an experimental processor is 
not expected to be identical to the Testbed solutions due to potentially different amounts 
of round-off error caused by the following factors. 

1. The coefficient matrix of the linear system to be solved by an experimental processor 
is ordered differently. That is, processors INV and SSOL solve (in double precision) 

(PAP r ) Px = Pf , 

whereas our experimental processor solves (in double precision) 

(PAP T ) Px = Pf . 

Since the permutation matrix P is induced by resequencing the joints in the model, 
it is not the same as the permutation matrix P chosen by SPARSPAK-A for the 
coefficient matrix in general. 
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2. Even for the same ordering of A, the factorization algorithms implemented by different 
processors employ a different computation sequence. 

3. The system stiffness matrix is ill-conditioned. 

However, with the condition number estimates available for each system stiffness matrix, 
we can estimate the relative error in our solution x by the algorithm described in reference [1] 
and implemented in SPARSPAK-A. On the other hand, by assuming that the Testbed 
solution x is the correct solution we cam compute the relative error in x by 


Vi \£i\ 

We can now verify the correctness of our experimental processors if the relative error com- 
puted above is very close to the relative error estimated by SPARSPAK-A with respect to 
the true (but unknown) solution. We have listed these two quantities in Table 19 for all test 
problems and we see that they are essentially of the same magnitude or sufficiently close 
for all problems. 

Table 19: Comparing NICE/SPAR solutions x with SPARSPAK-A solutions £. 


problem 

max 

SPARSPAK-A estimate of 
the relative error in x 

max “i et 

p648 

5.9 x 10- 8 

6.9 x lO" 9 | 

focusl 

6.4 x 10" 8 

2.7 x 10~ 8 

focus 2 

5.9 x 10- 8 

4.6 X lO" 8 

focus3 

5.8 x 10" 8 

4.9 X lO" 8 

focus4 

7.3 x 10- 8 

3.8 x lO" 8 

demol 

4.8 X 10- 8 

1.4 x lO" 8 

demo2 

5.0 x 10- 8 

5.4 x lO" 9 

demo3 

4.7 x lO" 7 

2.9 x lO' 8 

demo4 

5.0 x lO" 8 

6.4 x lO" 9 

demo6 

1.6 x lO" 6 

1.2 x 10“ 8 

demo 7 

6.2 x lO" 8 

1.9 x 10~ 8 | 

demo9 

2.7 x 10- 8 

1.7 x 10~ 9 | 

demol 0 

5.6 x lO" 6 

1.6 x 10" 5 

demol 2 

5.7 X lO" 8 

4.4 x 10- 6 

demol 3 

5.8 x lO" 7 

6.9 x lO" 8 
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5*3 The Experimental Factorization Processors 

In this section, we briefly describe the three sparse matrix factorization processors we have 
installed in the CSM Testbed. The three processors employ different methods in solving a 
sparse symmetric positive definite system 


Ax — b . 

1. Processor SPK: The method employed by the processor SPK is the direct solver pro- 
vided by SPARSPAK-A corresponding to the minimum degree ordering algorithm in 
reference [19]. 

2. Processor EXP1: The factorization method employed by the experimental processor 
EXP1 is the multifrontal method implemented by Liu as described in reference [16]. 

3. Processor EXP2: The factorization method employed by the experimental processor 
EXP2 is the adaptive sparse out-of-core Cholesky scheme recently developed by Liu 
[15]. 

Since the factorization methods employed by the processors EXP1 and EXP2 use the same 
storage scheme as that used by the minimum degree ordering in SPARSPAK-A and they 
were intended to be used :in conjunction with SPARSPAK-A [15, 16], the same interface 
modules for inputing the problem to SPARSPAK-A can be used. 

5.4 Numerical Results 

We first compare the factorization time of the three experimental processors with that of 
the processor INV. Since the joint ordering can affect the execution time of processor INV 
significantly, we have attempted to report the timing results for all available joint elimination 
sequences. The ordering algorithms currently available in the CSM Testbed are listed in 
Table 20. 


Table 20: The joint ordering methods employed in the CSM Testbed. 


acronym 

ordering algorithm 

ND 

Nested dissection (fill minimizer) [12] 

MDG 

Minimum degree (fill minimizer) [12, 17] 

RCM 

Reverse Cuthill-McKee (profile minimizer) [12] 

GPS 

Gibbs-Poole-Stockmeyer (bandwidth minimizer) [3] 

SEQ 

Sequential joint elimination sequence (i.e., no reordering of joints) 


Since the ordering algorithms used by processors EXP1 and EXP2 are the topological 
orderings of the elimination tree induced by the minimum degree ordering [15, 18], we have 
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thus used “MDG*” to represent any one of them. One consequence of the choice of ordering 
algorithms by the experimental processors is that the amount of fill-in in the Cholesky factor 
is the same for the three of them. From the factorization times reported in Table 21 we 
see that the in-core factorization time of processors SPK and EXP1 are significantly smaller 
than the INV times in most cases as one would expect in view of the I/O conducted by 
the latter. Except for problem demo7, the saving in execution time ranges from 30% to 
58% compared to the fastest INV time. As we have pointed out earlier, the reordering of 


Table 21: NICE/SPAR (INV) and SPARSPAK-A factorization times. 


Factorization times (in seconds) 


NICE/SPAR (INV) 

SPK 

EXP1 

EXP 2 



ND 

MDG 

RCM 

GPS 

MDG* 

MDG* 


focusl 

466 



196 

239 

44 

43 

■ai 

focus 2 


236 

145 

445 


76 

76 

102 

focus3 


770 

441 



■a 


378 





940 


148 

146 

188 

demo6 

61 

53 

mi 

82 

82 

33 

32 

45 

demo 7 

96 

124 


113 

112 

93 

92 

mm 

demo 10 

60 

304 

62 

203 

166 

41 

40 

HK9I 

demo 13 

406 





283 




the joint sequence in the model produces a different permutation matrix from that induced 
by applying the same ordering algorithm to the coefficient matrix itself. In Table 22, we 
compare the quality of the minimum degree algorithm when applying to each case, where 
we give the nonzero counts in the system stiffness matrix A and the computed factors. Due 
to the different storage schemes employed by the processor INV and SPARSPAK-A, the 
fill-in is not measured in exactly the same manner as Table 22 indicates. 


Table 22: Comparing the fill-in of different processors. 



Processor INV 

SPK, EXP1, EXP2 ~j] 

problem 

ordering 

# nonzeros 
in K.SPAR 



# off-diag 
nonzeros in A 

# off-diag 
nonzeros in L 


MDG 


71040 


18458 

47487 I 


MDG 

45792 


MDG* 

26824 

72519 


MDG 



MDG* 


184042 

focus4 

MDG 

65088 

165480 

MDG* 

38656 

120682 

demo 7 

MDG 

22464 


MDG* 

18576 

62829 

demo 10 

MDG 

47376 

79992 

MDG* 


71076 
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The performance of processors SPK and EXP1 are essentially the same in terms of ex- 
ecution time. In terms of storage, the in-core multifrontal Cholesky factorization scheme 
of processor EXP1 requires additional working storage compared with the in-place Cholesky 
method provided by processor SPK. However, it should be pointed out that the multifrontal 
method lends itself readily to out-of-core implementation [20], in which case the amount 
of in-core storage required to perform the entire factorization turns out to be precisely the 
same as the required working storage for the in-core version. Readers are referred to ref- 
erence [20] for various strategies in minimizing the working storage. In reference [18] the 
behaviour of the multifrontal method in a paging environment is studied. 

In order to compare the out-of-core performance of processor EXP2 with that of processor 
INV, we should note the following. 

1. The number of in-core data reorganizations of the adaptive sparse out-of-core Cholesky 
scheme [15] is dynamically adjusted to the available memory. In particular, if the de- 
clared working space is sufficiently large for the given problem, the entire factorization 
process will be carried out in-core without reorganizing the data structures. In order 
to provide a meaningful comparison of the performance of processor EXP2 in execution 
time as well as storage requirement with that of processor INV, we have run the proces- 
sor EXP2 with the minimum amount of in-core storage that will allow EXP2 to execute. 
This number can be determined in advance of the actual numerical factorization. 

2. The processor EXP2 does I/O using ordinary text files. In particular, the sparse 
coefficient matrix is saved in a text file and read into memory one column at a time, 
and the computed Cholesky factor is written to a text file one column at a time. In the 
current implementation, auxiliary storage is not used to reduce the in-core overhead 
storage, although it is possible to do so as suggested in reference [15]. 

3. We have explained in detail how the processor INV carries out the out-of-core block 
LDL t factorization process in Section 3.2 of this report. The I/O traffic involved 
amounts to retrieving the system stiffness matrix from the dataset K.SPAR .* as well 
as the indexing information from the dataset AMAP..ic2.isize , and outputting the 
computed factors to the dataset INV.K.ncon. Because the data are read from or 
written to the database one record at a time, the number of disk I/O operations is 
determined by the record length of each dataset. The default record length of these 
three datasets are listed below in Table 23. 

Recall that one record has to accommodate at least the amount of data needed to pro- 
cess one block row of the coefficient matrix. Hence the default record length may not 
be big enough for larger or denser problems and they can again be changed by reset- 
ting the designated argument when executing the source processor of each respective 
dataset. In particular, if necessary, processor T0P0 will automatically increase the 
AMAP record length twice up to a maximum size of 2.25xLRAMAP words. The 
number of records contained in each dataset are given under the column heading 
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Table 23: Database interface of processor INV. 


Database interface of processor INV 

Source processor 

Reset argument 

Dataset name 

Default record length 

K 

LREC 

K.SPAR.* 

2240 words 

TOPO 

LRAMAP 

AMAP..ic2.isize 

1792 words 

INV 

LRA 

INV.K.ncon 

3584 words 


“Records” in the table of contents of the data library created for each particular 
analysis. 

In summary, the volume of I/O involving each individual dataset is roughly the prod- 
uct of the number of records and the record length (strictly speaking, the last record 
may contain fewer items than axe permitted by the specified record length), whereas 
the number of disk read/write operations is determined by the number of records. 

4. The in-core storage required by the processor INV must accommodate one record 
of each dataset in Table 23 in addition to accommodating the maximum number 
of submatrices involved during the factorization process. Therefore, as suggested in 
reference [24], the memory requirement for processor INV may be estimated by the 
following formula. 

number of words = J + X 3 + m ( L\ + L 2 + n 2 /*) , 


where 

J = the number of joints in the structure. 

L\ = record length of input dataset K.SPAR.jdf2. 

L 2 — record length of INV.K.ncon dataset. 

L$ — record length of AMAP..ic2.isize dataset. 
m = 1 for single precision; 2 for double precision. 

n = maximum number of degrees of freedom per joint (default 3, 4, 5, or 6). 

I s = the maximum number of submatrices in use during any one stage of the fac- 
torization process. Its value can be obtained from the processor TOPO output 
parameter SIZE INDEX or from the value of isize from A MA P. . ic2. isize. 

It was suggested in reference [24] that this formula may be used to estimate the 
amount of space in blank common required by processor INV. If the number of words 
required is larger than the dimension of blank common, the blank common dimension 
must be increased and the Testbed must be recompiled. 
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In Table 24, we compare the factorization time and the memory requirement of processor 
INV with that of the experimental processor EXP2. For each problem, we give the number 
of nonzero elements in the Cholesky factor computed by SPARSPAK-A (recall that the 
amount of fill-in is the same for all three experimental processors) under the column heading 
“NOFNZ”. The ratio of memory requirement to the size of the computed Cholesky factor 
is computed for each problem and displayed for both processors. Note that the quantity of 
n 2 I s we use in measuring the memory requirement of processor INV is an underestimate as 
explained above. We use “LNZSZE” to indicate the maximum number of nonzeros which 
have to be present in-core for the adaptive sparse Cholesky factorization process to be 
successfully executed. The results in Table 24 indicate that the processor EXP2 can be quite 
competitive in both time and space. 


Table 24: Comparing two out-of-core factorization processors. 




INV 

EXP2 

problem 

NOFNZ 

n 2 Ij 

Time (sec) 

LNZSZE 

Time (sec) 

NOFNZ 

7TUFWZ 

focusl 

47487 

61% 

107 

35% 

61 

focus2 

72519 

47% 

147 

31% 

102 

focus3 

184042 

36% 

449 

32% 

378 

focus4 

120682 

46% 

288 

31% 

188 

demo6 

28302 

62% 

53" 

44% 

45 

demo 7 

62829 

46% 

108 

45% 

113 

demolO 

71076 

39% 

62 

5% 

72 

demo 13 

180315 

6% 

406 

25% 

337 


Comparing the factorization algorithm of processor EXP2 (adaptive out-of-core Cholesky) 
with that of processor SPK, we see that the difference in their execution time can be ac- 
counted for in the following three aspects. 

1. The time spent in data structure reorganization. 

2. The time for reading in the coefficient matrix A column by column. 

3. The time for writing out the computed Cholesky factor L column by column. 

The timing results reported in Table 24 are those with the minimum amount of mem- 
ory and maximum number of data structure reorganizations. Since the frequency of data 
structure reorganizations can be reduced by providing more memory, there is a potential 
tradeoff between time and space. However, the timing results in Table 25 indicate that the 
time spent in this regard is too small to justify the more significant increase in storage. We 
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Table 25: Data structure reorganization and factorization time. 


n 

EXP 2 

EXP 2 

problem 

NOFNZ 

LNZSZE 

# REORGZ 

Time 

LNZSZE 

# REORGZ 

Time 

focus 1 

47487 

28872 

2 

59 sec 

16823 

29 

61 sec 

focus 2 

72519 

34020 

3 

99 sec 

22647 

31 

102 sec 

demo 6 

28302 

17676 

3 

44 sec 

12394 

34 

45 sec 

demo 10 

71076 

27900 

2 

69 sec 

3330 

169 

72 sec 


can thus conclude that the I/O time can be considered to be the sole factor in determining 
the speed of processor EXP2. 

Since the multifrontal Cholesky method is also a good candidate for out-of-core imple- 
mentation, and we pointed out earlier that the “working storage” required in its in-core ver- 
sion is precisely what is needed as working storage in its out-of-core version, it makes sense 
to evaluate its out-of-core potential by comparing its minimum working storage requirement 
with the memory requirement of processor EXP2. The results we present in Table 26 indicate 
that the two are quite comparable as far as the test problems are concerned. 


Table 26: Comparing processor EXP1 with EXP2 



Multifrontal (EXPl) 

Column- Cholesky (EXP2) 

problem 

NOFNZ 

LNZSZE/NOFNZ 

LNZSZE/NOFNZ 

focusl 

47487 

34% 

35% 

focus2 

72519 

31% 

31% 

focus3 

184042 

36% 

32% 

focus4 

120682 

28% 

31% 

demo6 

28302 

44% 

44% 

demo 7 

62829 

56% 

45% 

demolO 

71076 

6% 

5% 

demol3 

180315 

21% 

25% 


For completeness, we provide in Table 27 the timing results of three other processors 
which are also essential in solving the linear system arising from a Testbed problem, namely 
TOPO, K and SSOL. 

Finally, we provide in Table 28 the total time in executing the processor SPK in the 
Testbed and indicate separately the time attributed to the numerical factorization phase 
and the triangular solution phase. The SPK time thus includes the time for retrieving data 
from the global database and setting up the problem for the SPARSPAK-A solver. 

In summary, our preliminary findings indicate that there are alternative sparse matrix 
techniques which are suitable for more general applications and appear to be also competi- 
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Table 27: Timing results of TOP0 , K, INV, SSOL. 


problem 

TOPO 




EESSSSI 

focusl 

5 




148 sec 

focus2 

6 

34 


Ea 

204 sec 

focus3 



EH 


559 sec 

focus4 



288 

HI 

inn 

de:mo6 

3 


53 

9 

mam 


4 


ESI 


§|y0gg|j 

demolO 

5 

18 


16 

101 sec 

demol3 

11 


EH 


HI 


Table 28: Execution time of the processor SPK. 


problem 

fact 


h m 

IKBI 


2 sec 





El 

focus 3 



376 sec 

focus4 



194 sec 

demo6 

33 sec 

2 sec 

48 sec 

demo7 



H 

demolO 

EEEH 


73 sec 

demo 13 

283 sec 

9 sec 

331 sec 
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tive in execution time and storage usage compared to the techniques currently employed in 
the CSM Testbed. 
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Appendix A Installing the Processor SPK 

The processor SPK consists of a subset of SPARSPAK-A [2] modules and a set of subroutines 
which provide an interface between SPARSPAK-A and the global database of the CSM 
testbed. All of the subroutines are provided as a single directory SPARSE on a UNIX tax 
tape. The Fortran source for the package is distributed among a number of subdirectories. 

There are “make” files provided, so that the person installing the package needs only to 
execute a few commands to compile the package and create the run-time library. 

It is advisable to read “§4 Developing New Matrix Factorization Processors” of this 
report before beginning installation of the package. Since the SPARSE package is used in 
conjunction with the CSM testbed, we assume in the sequel that the NICE/SPAR processors 
have been properly installed in the directory /usr/ns/nice and /usr/ns/spar, and that 
the SPARSE package is to be installed in the directory /usr/ns/sparse. The hierarchy 
of the directory /usr/ns and the files relevant to the installation and use of the SPARSE 
package are depicted in Figure 30. 

The steps to install the SPARSE package are as follows. 

1. Create a directory for SPARSE: 

cd /usr/ns 
mkdir sparse 
cd sparse 

2. Copy the files from tape to disk: Put the tape in the tape drive and “tar” the files to 
the new disk directory: 

tar xvf /d e\/ device 

where device should be the appropriate name of the tape drive on your machine. Do 
an “Is” to make sure that three directories (install, csm-intrface and spk-subset) 
have been copied from the tape. 

3. Edit the installation- dependent subroutines : The package has installation- dependent 
subroutines SPK, CTIME, SPK CSM, DTIME and SPRSPK which provide timing 
information to the package and set some installation-dependent parameters. In ap- 
pendix §B, we provide a set of examples for these subroutines. The sample programs 
are written for a SUN/3 workstation running the UNIX operating system at the Uni- 
versity of Tennessee Knoxville. Comments in these listings indicate changes which may 
be necessary. The subroutine SPK is contained in the directory csm-intrf ace/driver, 

the subroutines CTIME and SPKCSM are contained in the directory csm-intrf ace/system, 
and DTIME and SPRSPK are contained in the directory spk-subset/system. Sam- 
ples of subroutines required by CTIME, DTIME and SPRSPK can be found in the 
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/usr/ns 


nice spar sparse 


install sparselib.a spk-subset csm-intrface 


Makefile nicespar .ams makefile .ns . spar driver 






spk.f spka.f spkobjs.a 


Figure 30: The file system of the directory /usr/ns. 
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directory spk-subset/local; these subroutines are appropriate for machines running 
Berkeley 4.2 or 4.3 UNIX and their derivatives such as SUN 05. 

4. Edit the make file /usr/ns/sparse/install/Makef ile : Compilation of the package 
is performed using a collection of UNIX make files. The most important make file is 
called Makefile found in the directory install; it will invoke the other make files. 
The distributed make files assume that the package is running on a SUN workstation. 
There are comments in Makefile to help you make the appropriate changes to it 
for your installation. There is no need to change the make files in any other sparse 
directories. 

5. Create and install the compiled library: After making the required changes to Makefile, 
you are ready to create and install the compiled library. Execute the following com- 
mands. 

cd /usr /ns /sparse/install 
make install 

A compiled library sparselib.a will be created in the directory sparse. 

6. Install a new processor in the testbed: Since the SPARSE package is installed as 
a processor SPK in the testbed and a CSM processor is a subroutine called by the 
NICE/SPAR main program, it is necessary to compile the SPK driver routines in the 
directory /usr/ns/sparse/csm-intrf ace/driver and edit the main program master 
file nicespar.ams in the directory /usr/ns/spar. The object code of the SPK driver 
routines spk.f and spka.f is contained in a separate library called spkobjs.a in 
the driver directory so that it may be updated independent of sparselib.a. In 
addition, the makefile in the directory /usr/ns/spar must be edited so that the two 
libraries can be linked to the executable when it is created. A copy of the properly 
edited nicespar.ams and a copy of the edited makefile can be found in the directory 
install. The former has the file name nicespar . ams and the latter has the file name 
makefile. ns. spar. With these two files available, the following commands may be 
executed to install the new processor SPK in the testbed. Note that you must have 
write permission in the directory spar to do this. 

cd /usr/ns/sparse/install 
make spk 
cd ../../spar 

mv makefile makefile. old 

mv nicespar.ams nicespar.ams.old 

cp ../sparse/install/makefile.ns.spar makefile 
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cp ../sparse/install/nicespar.ams 
make 


mcespar.ams 


7. When the file korcoma.inc is changed: Since the include file korcoma.inc in the 
directory spar declares the size of the in-core storage available for every SPAR pro- 
cessor, the driver source code spk.f of processor SPK must contain the line 

include ’/usr/ns/spar/korcoma.inc’ 

and it must be recompiled each time the declared size is changed. Since the depen- 
dence of spk . o on korcoma . inc is specified in the appropriate make file, the following 
commands will not only detect whether the declaration file korcoma.inc has been 
modified since spk.o was last created but also recompile spk.f and update the li- 
brary spkobjs.a if that is the case. Finally the executable in the directory spar is 
recreated to link to the modified spkobjs.a after the “make” command in the last 
line is executed. 

cd /usr/ns/sparse/install 
make spk 
cd ../../spar 
make 

8. Recover space used by intermediate files: If the system on which you are running 
is short of disk space, a substantial amount of space used during the installation of 
SPARSE can be recovered by deleting the w .o” files and other intermediate files gen- 
erated during the creation of the library. To do this, execute the following commands. 

cd /usr/ns/sparse/install 
make clean 

If for some reason you must later re-create some or all of the library sparselib.a, 
these intermediate files will have to be regenerated, at considerable cost in computer 
time. Thus, it is advisable to execute “make clean” only if you really need the space. 


76 



ooooo oooo o oooooooo o a o o o a o o o o o a a a o o o o o o o o o o o o o o o o o o a a o o a 


Appendix B Installation-dependent Subroutines 


£****************************************************** *m******** ************ 
£ **** ****** ****** ****** * **** * * *** * ***** * * * * *** * * ********* * ***** * ***** *** w* * *** * 

C SPK A NEW CSM PROCESSOR 

£5 **** *** ** * * * ** * * ******************** ***** ****** * ** ********* **** ***** ** * ****** * 
£<**** ********* ***************************** ******** ********* ******************* 


PURPOSE - THIS IS THE DRIVER FOR INSTALLING INTO NICE/SPAR 
OUR INTERFACE MODULES AS A SINGLE PROCESSOR WHICH 
SOLVES CSM TESTBED PROBLEMS USING SPARSPAK-A MODULES. 

THE NEW PROCESSOR IS CODED AND INSTALLED INTO NICE/SPAR DIRECTLY 
FOLLOWING THE GUIDELINES GIVEN IN NASA TECHNICAL MEMORANDOM 
89096, NAMELY 

(a) THE NAME OF THE PROCESSOR SHOULD BE NO LONGER THAN FOUR 
CHARACTERS. 

(b) THE PROCESSOR SHOULD BE WRITTEN AS A FORTRAN 77 SUBROUTINE 
WHOSE NAME IS THE PROCESSOR NAME. 

(c) THE SUBROUTINE SHOULD HAVE NO ARGUMENTS. 

(d) THE PROCESSOR SHOULD BEGIN EXECUTION WITH A CALL TO THE 
LIBRARY SUBROUTINE “INTRO” WITH THE PROCESSOR NAME 

AS THE ONLY ARGUMENT. THE GIVEN NAME IS USED BY THE 
“GAL” DATA MANAGER AS THE CREATING PROCESSOR FOR 
NEW DATASETS INSERTED IN “GAL” LIBRARIES; IT ALSO 
APPEARS IN THE INTERACTIVE PROMPT STRING IF THE 
“SPAR READER” ROUTINE IS USED FOR INPUT COMMAND 
PROCESSING. 

(e) THE LABELED COMMON BLOCK /I ANDO/ WITH 2 INTEGER VARIABLES 
CONTAINING USER INPUT AND OUTPUT UNIT NUMBERS SHOULD BE 
INCLUDED IN APPROPRIATE MODULES. THE UNIT NUMBERS ARB 
ASSIGNED IN THE SUBROUTINE “INTRO”. 

(f) CALL LIBRARY SUBROUTINE “FIN” TO CLOSE “GAL” LIBRARIES. 

*************************************************************** 

WARNING 

*************************************************************** 

THE PATHNAME OF THE INCLUDE FILE “korcoma.iac” 

IS INSTALLATION DEPENDENT. 


***************************************************************************** 


SUBROUTINE SPK 


INCLUDE DECLARATION CONTAINING BLANK COMMON VARIABLES AND 
DIMENSIONS: 

PARAMETER (KSZZZ= 200000) 

COMMON KORB, KEVBN, KORT, A(KSZZZ) 


include ’/ n,r *MC68020/ nlal } echu/m/ tpar/ korcoma.inc’ 
INTEGER MX STORE 


IDENTIFY PROCESSOR TO CSM ARCHITECTURE 
CALL INTRO ( ’SPK’ ) 

WORKING STORAGE A IS DECLARED AS KSZZZ WORDS WHICH IS 
EQUIVALENT TO HALF THAT MANY DOUBLB-PRECISION FLOATING 
POINT NUMBERS. 


MXSTOR = KSZZZ/2 
CALL SPKA ( A, MXSTOR ) 
CALL FIN (0,0) 

CALL EXIT 
END 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


0 ***************************************** ************************************* 
Q ***** **** ************* ***** ** mm* ********* *** **** *********** **** ********* ****** 

C CTIME ELAPSED PROCESSOR TIME 

Q * * *** * **** ****** ** **** **** * *** * * * ** * ** * 41 * * * * * ***** ***** m *** * Dc *** * ************ * 
0 ****************************************************************************** 
c 

C PURPOSE - CTIME RETURNS THE ELAPSED PROCESSOR TIME SINCE 
C IT WAS LAST CALLED. IT USES THE COMMON VARIABLE TIME 

C TO REMEMBER THE TIME WHEN CTIME WAS LAST CALLED. 

C 

0 ************************************************************************* 

C WARNING 

0 ************************************************************************* 

C THIS IS AN INSTALLATION DEPENDENT ROUTINE. IT 
C SHOULD BE SET UP BY THE INSTALLER OF THE PACKAGE. 

C IN THIS EXAMPLE, ROUTINE G TIMER IS THE TIMER ROUTINE 

C THAT RETURNS THE CURRENT PROCESSOR TIME ON A SUN/3 

C WORKSTATION RUNNING THE UNIX OPERATING SYSTEM AT THE 

C UNIVERSITY OF TENNESSEE KNOXVILLE. 

0 ************************************************************************* 

C 

C INPUT PARAMETER - 

C IDUMMY - A DUMMY INTEGER VARIABLE, 

C 

C PROGRAM SUBROUTINE - 

C GTIMER. 

C 

0 ************* ************************************* **«* ************************ 

c 

REAL FUNCTION CTIME ( IDUMMY ) 


<**************************************************************************** 

INTEGER IDUMMY, IPRNTE, IPRNTS, MAXINT 
REAL RATIOL, RATIOS, TIME ,X 

:****************************************************************** ********* * 


COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 

1 TIME 

****************************************************************************** 


CALL GTIMER ( X ) 
CTIME = X - TIME 
TIME = X 
RETURN 
C 

END 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


C******************************** ********************************************** 
£****************************************************************************** 
C SPKCSM INITIALIZE PARAMETERS 

£■**** ********* ******************* **** ********* ********* ***** ************* *41*41** 
£******************************** ********************** ************** ********** 


PURPOSE - TO SET SYSTEM PARAMETERS AND ASSIGN DEFAULT 

VALUES TO SOME USER PARAMETERS. IT IS A MACHINE 

DEPENDENT ROUTINE. THIS ROUTINE HAS TO BE CALLED 

BEFORE ANY OTHER PACKAGE MODULE. 

PARAMETERS INITIALIZED - 

IPRNTB . THE OUTPUT UNIT NUMBER FOR ERROR MESSAGES. 

IPRNTS - THE OUTPUT UNIT NUMBER FOR STATISTICS. 

RATIOL - THE RATIO OF THE NUMBER OF BITS IN A FLOATING 
POINT VARIABLE TO THAT IN A LONG INTEGER 
VARIABLE. FOR EXAMPLE, IF FLOATING POINT 
NUMBERS OCCUPY TWICE AS MANY BITS AS LONG 
INTEGERS, RATIOL SHOULD BE SET TO 2. 

RATIOS - THE RATIO OF THE NUMBER OF BITS IN A FLOATING 
POINT VARIABLE TO THAT IN A SHORT INTEGER 
VARIABLE. 

MAXINT - THE LARGEST POSITIVE INTEGER THAT CAN BE 
STORED IN A SHORT INTEGER VARIABLE. 

TIME - VARIABLE USED BY THE TIMER ROUTINE CTIME. 

SEE REMARK 

STAGE - STARTING STAGE OF SYSTEM-CSM. 

REMARK - THIS INTERFACE PACKAGE ASSUMES THE EXISTENCE OF 
A REAL TIME FUNCTION CTIME WHICH RETURNS THE ELAPSED 
PROCESSOR TIME SINCE IT WAS LAST CALLED. WITH THE 
COMMON VARIABLE TIME, THE INSTALLER OF THE PACKAGE 
SHOULD BE ABLE TO WRITE SUCH A FUNCTION, USING THE 
INSTALLATION TIMER. 


***** c***************** ********** ************* ***** ********* ******************* 


SUBROUTINE SPKCSM 

****************************************************************************** 


CHARACTER*^ LIBNAM 

CHARACTERS CDUMMY 

INTEGERM IIN, IOUTX 

INTEGERM IPRNTB, IPRNTS, MAXINT 

INTEGER** BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER** MSGLVL, IERR , MAXCSM 

REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTB, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 
COMMON /IANDO/ IIN, IOUTX 

*********************************************************** ************** 
WARNING 


THE FOLLOWING 4 LINES OF CODE ARE INSTALLATION 
DEPENDENT. THEY MAY HAVE TO BE MODIFIED BY THE 
PERSON INSTALLING THIS PACKAGE. 

OUR CURRENT ENVIRONMENT - 

. RATIOL AND RATIOS ARE BOTH 2. 

- MAXINT = 2**15 - 1 = 32767 


C 

C 

c 


INSTALLATION DEPENDENT PARAMETERS 

TIME = 0.0 

RATIOL s 2.0 
RATIOS = 2.0 

MAXINT = 32767 
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IPRNTE AND IPRNTS ARE BOTH SET TO THE WRITER UNIT 
NUMBER ASSIGNED TO IOUTX WHEN THE NEW PROCESSOR 
IS IDENTIFIED TO THE C SM- ARCHITECTURE. 


IPRNTE a IOUTX 
IPRNTS = IOUTX 


INITIALIZING THE EXECUTION STAGE FOR THE INTERFACE ... 


C 

c 


STAGE = 0 
RETURN 
END 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


C— SPARSPAK-A (ANSI FORTRAN) RELEASE III — NAME = DTIME 


C (C) UNIVERSITY OF WATERLOO JANUARY 1984 

Q ************* ***** ********* *************************************************** 
0 *»** ********* *********************** ********* ********* ************************ 
C DTIME DELTA TIME 

0 ********* ********* **** ************** ***** ****** ******* ****************** ****** 
0 ********* ********* *************************** *********************** **** ***** * 
C 

C PURPOSE . DTIME RETURNS THE ELAPSED PROCESSOR TIME SINCE 
C IT WAS LAST CALLED. IT USES THE COMMON VARIABLE TIME 

C TO REMEMBER THE TIME WHEN DTIME WAS LAST CALLED. 

C 

0 ************************************************************************ 

C WARNING 

0 *************************************************************** ********* 

C THIS IS AN INSTALLATION DEPENDENT ROUTINE. IT 

C SHOULD BE SET UP BY THE INSTALLER OF THE PACKAGE. 

C IN THIS EXAMPLE, ROUTINE GTIMER IS THE TIMER ROUTINE 

C THAT RETURNS THE CURRENT PROCESSOR TIME ON A SUN/3 

C WORKSTATION RUNNING THE UNIX OPERATING SYSTEM AT THE 

C UNIVERSITY OF TENNESSEE KNOXVILLE. 

0 ************************************************************************ 

C 

C INPUT PARAMETER - 

C IDUMMY - A DUMMY INTEGER VARIABLE. 

C 

C PROGRAM SUBROUTINE . 

C GTIMER. 

C 

0#««**I««*»«»*»*«**«I|I******************«*«***«» ********************************* 

C 

REAL FUNCTION DTIME ( IDUMMY ) 


umi**************************************************************************** 

INTEGER IDUMMY, IPRNTE, IPRNTS, MAXINT 
REAL MCHEPS, RATIOL, RATIOS, TIME , X 

****************************************************************************** 

COMMON /SPKSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 

1 MCHEPS, TIME 

****************************************************************************** 

CALL GTIMER ( X ) 

DTIME = X - TIME 
TIME = X 
RETURN 
C 

END 
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— SPARSPAK-A (ANSI FORTRAN) RELEASE III — NAME = SPRSPK 
(C) UNIVERSITY OF WATERLOO JANUARY 1984 
***ni«*ii(**************************ai******************** ************** ********** 
****************************************************************************** 
SPRSPK START SPARSPAK-A 

****************************************************************************** 

****************************************************************************** 

PURPOSE - TO SET SYSTEM PARAMETERS AND ASSIGN DEFAULT 
VALUES TO SOME USER PARAMETERS. IT IS A MACHINE 
DEPENDENT ROUTINE. THIS ROUTINE HAS TO BE CALLED 
BEFORE ANY OTHER PACKAGE MODULE. 

PARAMETERS INITIALIZED - 

IPRNTE - THE OUTPUT UNIT NUMBER FOR ERROR MESSAGES. 

IPRNTS - THE OUTPUT UNIT NUMBER FOR STATISTICS. 

RATIOL - THE RATIO OF THE NUMBER OF BITS IN A FLOATING 
POINT VARIABLE TO THAT IN A LONG INTEGER 
VARIABLE. FOR EXAMPLE, IF FLOATING POINT 
NUMBERS OCCUPY TWICE AS MANY BITS AS LONG 
INTEGERS, RATIOL SHOULD BE SET TO 2. 

RATIOS - THE RATIO OF THE NUMBER OF BITS IN A FLOATING 
POINT VARIABLE TO THAT IN A SHORT INTEGER 
VARIABLE. 

MAXINT - THE LARGEST POSITIVE INTEGER THAT CAN BE 
STORED IN A SHORT INTEGER VARIABLE. 

MCHEPS - THE MACHINE EPSILON (UNIT ROUNDOFF ERROR). 

TIME - VARIABLE USED BY THE TIMER ROUTINE DTIME. 

SEE REMARK. 

STAGEA - STAGE VARIABLE FOR SYSTEM-A. 

REMARK - THIS PACKAGE ASSUMES THE EXISTENCE OF A REAL TIME 
FUNCTION DTIME WHICH RETURNS THE ELAPSED PROCESSOR TIME 
SINCE IT WAS LAST CALLED. WITH THE COMMON VARIABLE 
TIME, THE INSTALLER OF THE PACKAGE SHOULD BE ABLE TO 
WRITE SUCH A FUNCTION, USING THE INSTALLATION TIMER. 

PROGRAM SUBROUTINES - 
ALLOW , STIMBR. 

****************************************************** ************************ 
SUBROUTINE SPRSPK 


INTEGER ICPADA, ICPADB, IERRA , IBRRB , IPRNTE, 

1 IPRNTS, MAXINT, MAXSA , MAXSB , MCOLS , 

1 MDCONS, MDEQNS, MSCONS, MSEQNS, MSGLVA, 

1 MSGLVB, NVARS , STAGEA, STAGEB 

INTEGER IIN, IOUTX 

REAL MCHEPS, RATIOL, RATIOS, TIME 
DOUBLE PRECISION EPS , EPS1 

:******************** *********************** ***** **** ***** **** ********* ***** * 


COMMON /SPKSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 
1 MCHEPS, TIME 

COMMON /SPAUSR/ MSGLVA, IERRA , MAXSA , NVARS 
COMMON /SPACON/ STAGEA, ICPADA(49) 

COMMON /SPBUSR/ MSGLVB, IBRRB , MAXSB , MCOLS , MSEQNS, 
1 MDEQNS, MSCONS, MDCONS 

COMMON /SPBCON/ STAGEB, ICPADB(49) 

COMMON /IANDO/ IIN, IOUTX 
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*** ************** ********** ********* ********* ************** **** ***** ********** 


****************** ********* ********* ***** ********* ************* ******* 
WARNING 

********************************************************************** 

THE FOLLOWING 6 LINES OF CODE ARB INSTALLATION 
DEPENDENT. THEY MAY HAVE TO BE MODIFIED BY THE 
PERSON INSTALLING THIS PACKAGE. 

ON A SUN/3 WORKSTATION AT THE UNIVERSITY OF TENNESEE KNOXVILLE: 
. STIMER IS THE ROUTINE TO START THE TIMER, 

- ALLOW IS THE ROUTINE TO ALLOW FOR A NUMBER OF 
ARITHMETIC UNDERFLOWS BEFORE SYSTEM ABORTS. 

- RATIOL AND RATIOS ARE 2 AND 4 RESPECTIVELY. 

. MAXINT = 2**15 - 1 = 32767 


TIME = 0.0 

CALL STIMER 

CALL ALLOW ( 1234567 ) 

RATIOL =2.0 
RATIOS = 2.0 

MAXINT = 32 767 


IPRNTE AND IPRNTS ARB BOTH SET TO THE WRITER UNIT 
NUMBER ASSIGNED TO IOUTX WHEN THE NEW PROCESSOR 
IS IDENTIFIED TO THE CSM- ARCHITECTURE. 


IPRNTE = IOUTX 
IPRNTS = IOUTX 


COMPUTE THE MACHINE EPSILON. 


BPS = 1.0D0 
LOO CONTINUE 

EPS = BPS/2. 0D0 
EPS1 = 1.0D0 + BPS 
IF ( BP Si .GT. 1.0D0 ) GO TO 100 
MCHBPS = BPS*2.0D0 

WRITE (IPRNTS, 11) 

11 FORMAT ( IH1 

1 /5X, 40H *** UNIVERSITY OF WATERLOO 

1 /5X, 40H********** SPARSE MATRIX PACKAGE 

1 /6X, 4 0H********** (SPARSPAK) 

1 / 5X, 40H* ********* RELEASE 3 

1 /5X, 4 OH********** (C) JANUARY 1984 ) 

WRITE (IPRNTS, 22) 

22 FORMAT ( 5X, ANSI FORTRAN ) 

WRITE (IPRNTS, 33) 

33 FORMAT ( 5X, 40 H********* * DOUBLE PRECISION ) 

WRITE (IPRNTS, 44) 

44 FORMAT ( 5X, 40 H********* * LAST UPDATE JANUARY 1984 ) 

WRITE (IP RNTS.55) IPRNTE, IPRNTS 
55 FORMAT (// 10X, 35HOUTPUT UNIT FOR ERROR MESSAGES , 17 
1 /10X, 35HOUTPUT UNIT FOR STATISTICS , 17 ) 


INITIALIZING USER VARIABLES FOR S YSTBM-A ... 


STAGEA = 0 
C 

RETURN 

C 

END 
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Appendix C Listing of Programs 


************* ********* ** *** ********* ********* *********** ********** *** ********* 
****#*iK**»«***)|i**«*W*«*»*«*****«»*»»****«************it»|i**********i|»l»l<**iK**»**** 

SPK A NEW CSM PROCESSOR 

********* ********* **** turn *** ***** **** **«*« * ******** **** ***** **** * **** **** **** * * 
*******************************************»*»*********************(*********** 


PURPOSE - THIS IS THE DRIVER FOR INSTALLING INTO NICE/ SPAR 
OUR INTERFACE MODULES AS A SINGLE PROCESSOR WHICH 
SOLVES CSM TESTBED PROBLEMS USING SPARSPAK-A MODULES. 

THE NEW PROCESSOR IS CODED AND INSTALLED INTO NICE/SPAR DIRECTLY 

FOLLOWING THE GUIDELINES GIVEN IN NASA TECHNICAL MEMORANDOM 

89096, NAMELY 

(*) THE NAME OF THE PROCESSOR SHOULD BE NO LONGER THAN FOUR 
CHARACTERS. 

(b) THE PROCESSOR SHOULD BE WRITTEN AS A FORTRAN 77 SUBROUTINE 
WHOSE NAME IS THE PROCESSOR NAME. 

(c) THE SUBROUTINE SHOULD HAVE NO ARGUMENTS. 

(d) THE PROCESSOR SHOULD BEGIN EXECUTION WITH A CALL TO THE 
LIBRARY SUBROUTINE “INTRO” WITH THE PROCESSOR NAME 

AS THE ONLY ARGUMENT. THE GIVEN NAME IS USED BY THE 
“GAL” DATA MANAGER AS THE CREATING PROCESSOR FOR 
NEW DATASETS INSERTED IN “GAL” LIBRARIES; IT ALSO 
APPEARS IN THE INTERACTIVE PROMPT STRING IF THE 
“SPAR READER” ROUTINE IS USED FOR INPUT COMMAND 
PROCESSING. 

(e) THE LABELED COMMON BLOCK /IANDO/ WITH 2 INTEGER VARIABLES 
CONTAINING USER INPUT AND OUTPUT UNIT NUMBERS SHOULD BE 
INCLUDED IN APPROPRIATE MODULES. THE UNIT NUMBERS ARE 
ASSIGNED IN THE SUBROUTINE “INTRO”. 

(f) CALL LIBRARY SUBROUTINE “FIN” TO CLOSE “GAL” LIBRARIES. 


WARNING 


THE PATH NAME OF THE INCLUDE FILE “korcome.inc” 
IS INSTALLATION DEPENDENT. 


********************** ********** ********* *************************** ********* 
SUBROUTINE SPK 


INCLUDE DECLARATION CONTAINING BLANK COMMON VARIABLES AND 
DIMENSIONS: 

PARAMETER (KSZZZ= 200000) 

COMMON KORB, KEVEN, KORT, A(KSZZZ) 


include */ UST -MC68020/nl»l/echu/n>/>p&r/korcom*.inc’ 
INTEGER MXSTORE 


IDENTIFY PROCESSOR TO CSM ARCHITECTURE 


CALL INTRO ( ’SPK' ) 


WORKING STORAGE A IS DECLARED AS KSZZZ WORDS WHICH IS 
EQUIVALENT TO HALF THAT MANY DOUBLE-PRECISION FLOATING 
POINT NUMBERS. 


MXSTOR = KSZZZ/2 
CALL SPKA ( A, MXSTOR ) 
CALL FIN (0,0) 

CALL EXIT 
END 
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£****************************************************************************** 
Q * * * * * * * ** * * ****** * * * ** * * * **** **** 111 *** * * ** ****** * ******* * ** * * ******** * * *** **** * 


C SPKA A DRIVER FOR INTERFACE MODULES AND SPARSPAK-A 

C ************* ******* <********************************************************** 
C ** * ****** ***** *** * ****************** ***** ******** * ******************* * ** * **** * 

c 

C PURPOSE - THIS IS THE DRIVER CALLING INTERFACE MODULES TO 
C SOLVE CSM TESTBED PROBLEMS USING SPARSPAK-A MODULES. 

C 

C INPUT PARAMETERS - 

C A - AN ARRAY OF MXSTOR DOUBLE-PRECISION FLOATING POINT 

C NUMBERS. 

C MXSTOR - SIZE OF ARRAY A IN DOUBLE-PRECISION FLOATING-POINT 

C NUMBERS. 

C 

C USER INPUT - 

C MSGLVL - MESSAGE LEVEL FOR INTERFACE MODULES. 

C MSGLVA - MESSAGE LEVEL FOR SPARSPAK-A MODULES. 

C BUFMAX - MAXIMUM BUFFER LENGTH ANTICIPATED. 

C LIB NAM - NAME OF THE DATA LIBRARY. 

C JDFSET - NAME OF DATASET JDF1.BTAB.1.8 

C KMAP - NAME OF DATASET KMAP.O.mnbs.kme 

C KSPAR - NAME OF DATASET K.SPAR.jdf2.0 

C CON - NAME OF DATASET CON.O.ncon.O 

C APPLF - NAME OF DATASET APPL.FORC.Uet.l 

C APPLM - NAME OF DATASET APPL.MOTI.Uet.l 

C STATD - NAME OF DATASET STAT.DISP.i*ct.ncon 

C 

C INTERFACE MODULES - 

C SPKCSM, LIBOPN, CTIME, SPACE , GETJDF, GETDOF, GTZERO , GTCOND, 

C GTMOTI, GBTIJ , GTFORC, GTNUM8, STATCS, GETSOL. 

C 

C SPARSPAK-A INTERFACE MODULES - 
C SPRSPK, ORDRB5 , SOLVES, EREST5, STATSA. 

C 

C LOGICAL READER UNIT NUMBER FOR USER INPUT - 41 
C 


£************* ***** ********* ********* *** *************** ********* ********* ****** 
c 

SUBROUTINE SPKA ( A, MXSTOR ) 

C 

DOUBLE PRECISION A(l) 

INTEGER MXSTOR 


CHARACTER* 40 LIBNAM 

CHARACTERS JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

INTEGER*4 IPRNTE, IPRNTS, MAXINT 

INTEGER*4 MSGLVL , IERR , MAXCSM 

INTEGER*4 DOF, BUF, MASK, KC, ICLQ, FCON, SPK 

INTEGER*4 BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER*4 MAXDOF, NEQNS, NUMJNT 

INTEGERS MSGLVA, IERRA , MAXSA , NVARS 

REAL GZTIME, G CTIME, GIJTIM, GFTIME, GMTIME.GNTIME, 

1 CSMTIM, CSMSTR 

REAL RATIOS, RATIOL, TIME 

INTEGER*4 SPACE 

REAL CTIME 

COMMON f CSMSYS / IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /CSMMAP / DOF, BUF, MASK, KC, ICLQ, FCON, SPK 
COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMDTA/ GZTIME, G CTIME, GIJTIM, GFTIME, GMTIME.GNTIME, 
1 CSMTIM, CSMSTR 

COMMON /PRBLEM/ MAXDOF, NEQNS , NUMJNT 

COMMON /SPAUSR/ MSGLVA, IERRA , MAXSA , NVARS 

e********************* ************** ****************** ************** **** ***** ** 

INTEGERS JLONG, NLONG, CSIZB 

INTEGERS IDUMMY, INDATA 

REAL RN, RNJNT, ROFFS, ROFFL, DUMMY 
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DOUBLE PRECISION RELERR, RELRES 


INITIALIZE SPARSPAK-A AND SYSTEM TIMER 


CALL SPRSPK 

INITIALIZE THE CSM.SPARSPAK INTERFACE PACKAGE 


CALL SPKCSM 
SET MSGLVL AS DESIRED 


INDATA = 41 

READ ( INDATA, 12 ) MSGLVL 
12 FORMAT (14) 


SET MSGLVA AS DESIRED 


READ ( INDATA, 12 ) MSGLVA 


SET MAXIMUM BUFFER LENGTH 


READ ( INDATA, 12 ) BUFMAX 


INPUT NAME OF LIBRARY AND DATASETS FOR GIVEN PROBLEM 


READ ( INDATA, 22 ) LIBNAM 
FORMAT( A40 ) 

READ ( INDATA, 32 ) JDFSET 
READ ( INDATA, 32 ) KMAP 
READ ( INDATA, 32 ) KSPAR 
READ ( INDATA, 32 ) CON 
READ ( INDATA, 32 ) APPLF 
READ ( INDATA, 32 ) APPLM 
READ ( INDATA, 32 ) STATD 
FORMAT( A51 ) 

OPEN THE LIBRARY 


CALL LIBOPN 


INITIALIZE THE TIMER 


DUMMY = CTIME(O) 
MXRBQD = BUFMAX 


SIZE OF STORAGE ARRAY 


MAX CSM = MXSTOR 


CHECK MAXCSM AGAINST MXRBQD 


IF ( SPACE ( IDUMMY ) .NE. 0 ) GO TO 9999 


RETRIEVE TOTAL NUMBER OF JOINTS AND STORE IN NJMJNT 


CALL GETJDF ( A ) 

COMPUTE FURTHER STORAGE REQUIREMENT 


ROFFS = RATIOS - 0.01 
ROFFL = RATIOL - 0.01 
RNJNT = NUMJNT + 1 

JLONG = IFIX((RNJNT+ ROFFL) /RATIOL) 

MXREQD = JLONG + BUFMAX 

IF ( SPACE ( IDUMMY ) .NE. 0 ) GO TO 9999 


COMPUTE ADDRESSES 


DOF = 1 

BUF = DOF + JLONG 


RETRIEVE DEGREES OF FREEDOM PER JOINT, 
AND INITIALIZE MAXDOF AND NEQNS 


CALL GETDOF ( A(DOF), A(BUF) ) 
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ADJUST BUFFER SPACE 

MXREQD = MXREQD . BUFMAX 
BUFMAX = MAXO ( BUFMAX, NEQNS ) 

MXREQD = MXREQD + BUFMAX 

IF ( SPACE ( IDUMMY ) .NE. 0 ) GO TO 9999 

MXUSED s MXREQD 

COMPUTE FURTHER STORAGE REQUIREMENT 


RN = NEQNS 

NLONG = IFIX ((RN + ROFFL)/RATIOL) 

MXREQD = MXUSED + NLONG 

IF ( SPACE ( IDUMMY ) .NE. 0 ) GO TO 9999 

COMPUTE ADDRESSES 

MASK = BUF + BUFMAX 

DETECT DUMMY ROWS 

CALL GTZERO ( A(DOF), A(BUF), A(MASK) ) 

MXUSED = MXREQD 

COMPUTE FURTHER STORAGE REQUIREMENT 
MXREQD = MXUSED + 7 

IF ( SPACE ( IDUMMY ) .NE. 0 ) GO TO 9999 
COMPUTE ADDRESSES 
KC = MASK + NLONG 
DETECT CONSTRAINED VARIABLES 

CALL GTCOND ( A(DOF), A(BTJF), A(KC), A(MASK), CSIZE) 

MXUSED = MXREQD 

COMPUTE FURTHER STORAGE REQUIREMENT 

MXREQD = MXUSED + MAXDOF + CSIZE 
IF ( SPACE ( IDUMMY ) .NE. 0 ) GO TO 9999 

TOTAL STORAGE TO BE USED 

MXUSED = MXREQD 

COMPUTE ADDRESSES 

ICLQ =KC + 7 

FCON = ICLQ + MAXDOF 

GATHER NONZERO CONSTRAINTS 

CALL GTMOTI ( A(BUF), A(MASK), A(FCON), CSIZE ) 

INTERFACE WITH SPARSPAK-A 

SPK = MXUSED + 1 
MAXSA = MAXCSM • MXUSED 

INPUT NONZERO STRUCTURE TO SPARSPAK-A 

CALL GETIJ( A(DOF), A(BUF), A(ICLQ), A(MASK), A(SPK)) 

DETERMINE SYMMETRIC ORDERING 

CALL ORDRB5 ( A(SPK) ) 

INPUT RIGHT HAND SIDE 

CALL GTFORC( A(BUF),A(MASK), A(SPK) ) 

INPUT MATRIX COEFFICIENTS AND RIGHT HAND SIDE MODIFICATIONS 
CALL GTNUM5( A(DOF), A(BUF), A(MASK), A(FCON), A(SPK)) 
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PERFORM NUMERICAL FACTORIZATION AND SOLUTION 


CALL SOLVE5 ( A(SPK) ) 

CSMTIM = CTIME(O) 

CALL EREST5 ( RELERR, A(SPK)) 


COMPARE WITH KNOWN NICESPAR SOLUTION 


CALL GETSOL ( A(BUF), A(SPK), RELRES ) 
CALL STATCS 
CALL STATSA 
C 

9999 CONTINUE 
RETURN 
C 

END 
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Q **#* ***** **** ********* ***** ************** **** ******* dc m ***** * # ** * **** **** *m *** * 
Q ***** i************** ******************* ****************** **** ********* ****** 

C SPKCSM INITIALIZE PARAMETERS 

0****************** ********* ****** *************************************** ##**** 
Q#*###******** *###*#*************#***#*#***##*##*****##*****####***#***####***# 

c 

C PURPOSE - TO SET SYSTEM PARAMETERS AND ASSIGN DEFAULT 
C VALUES TO SOME USER PARAMETERS. IT IS A MACHINE 

C DEPENDENT ROUTINE. THIS ROUTINE HAS TO BE CALLED 

C BEFORE ANY OTHER PACKAGE MODULE. 

C 

C PARAMETERS INITIALIZED - 

C IPRNTE - THE OUTPUT UNIT NUMBER FOR ERROR MESSAGES. 

C IPRNTS - THE OUTPUT UNIT NUMBER FOR STATISTICS. 

C RATIOL - THE RATIO OF THE NUMBER OF BITS IN A FLOATING 

C POINT VARIABLE TO THAT IN A LONG INTEGER 

C VARIABLE. FOR EXAMPLE, IF FLOATING POINT 

C NUMBERS OCCUPY TWICE AS MANY BITS AS LONG 

C INTEGERS, RATIOL SHOULD BE SET TO 2. 

C RATIOS - THE RATIO OF THE NUMBER OF BITS IN A FLOATING 

C POINT VARIABLE TO THAT IN A SHORT INTEGER 

C VARIABLE. 

C MAXINT - THE LARGEST POSITIVE INTEGER THAT CAN BE 

C STORED IN A SHORT INTEGER VARIABLE. 

C TIME - VARIABLE USED BY THE TIMER ROUTINE CTIME. 

C SEE REMARK 

C STAGE - STARTING STAGE OF SYSTEM-CSM. 

C 

C REMARK - THIS INTERFACE PACKAGE ASSUMES THE EXISTENCE OF 
C A REAL TIME FUNCTION CTIME WHICH RETURNS THE ELAPSED 

C PROCESSOR TIME SINCE IT WAS LAST CALLED. WITH THE 

C COMMON VARIABLE TIME, THE INSTALLER OF THE PACKAGE 

C SHOULD BE ABLE TO WRITE SUCH A FUNCTION, USING THE 

C INSTALLATION TIMER. 

C 

C 

0 **# * * * * * *# * * * ****** * * ** ** ************ * ** * * * * **#*##***#****# * ***** * * ***** * **** * 

c 

SUBROUTINE SPKCSM 
C 

0 ********************** ************** ********* ***** ********* **** **##*#*** ***** * 

c 

CHARACTER* 40 LIBNAM 
CHARACTER* 51 CDUMMY 
INTEGER*4 IIN, IOUTX 
INTEGER*! IPRNTE, IPRNTS, MAXINT 
INTEGER*! BUFMAX, MXUSED, MXREQD, STAGE 
INTEGER*! MSGLVL, IERR , MAXCSM 
REAL RATIOS, RATIOL, TIME 
C 

COMMON / CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON / CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(T) 
COMMON /IAND O / IIN, IOUTX 
C 

q ******************************** **** * ******** 1*1 **** ** *** * * ***** ****** * **** 

C WARNING 

0 ************************************************************************* 

C 

C THE FOLLOWING 4 LINES OF CODE ARE INSTALLATION 

C DEPENDENT. THEY MAY HAVE TO BE MODIFIED BY THE 

C PERSON INSTALLING THIS PACKAGE. 

C 

C OUR CURRENT ENVIRONMENT - 

C - RATIOL AND RATIOS ARE BOTH 2. 

C - MAXINT = 2**15 . 1 = 32767 

C 

C 

C INSTALLATION DEPENDENT PARAMETERS 

C 

TIME = 0.0 
C 

RATIOL = 2.0 
RATIOS = 2.0 
C 

MAXINT = 32 767 
C 
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IPRNTE AND IPRNTS ARE BOTH SET TO THE WRITER UNIT 
NUMBER ASSIGNED TO IOUTX WHEN THE NEW PROCESSOR 
IS IDENTIFIED TO THE CSM-ARCHITECTURE. 


IPRNTB = IOUTX 
IPRNTS s IOUTX 


INITIALIZING THE EXECUTION STAGE FOR THE INTERFACE 


C 

C 


STAGE = 0 
RETURN 
END 
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Q **** ********** **** ********* ***** ********* ********* *************************** * 
0 ****************** ************************************************************ 


C CTIMB ELAPSED PROCESSOR TIME 

0 ******************************************************************** ********** 
^ ****************************************************** ************************ 
C 

C PURPOSE - CTIMB RETURNS THE ELAPSED PROCESSOR TIME SINCE 
C IT WAS LAST CALLED. IT USES THE COMMON VARIABLE TIME 
C TO REMEMBER THE TIME WHEN CTIMB WAS LAST CALLED. 

C 

q ****************** ************** ***************************************** 

C WARNING 

0 ************************************************************************* 

C THIS IS AN INSTALLATION DEPENDENT ROUTINE. IT 

C SHOULD BE SET UP BY THE INSTALLER OF THE PACKAGE. 

C IN THIS EXAMPLE, ROUTINE G TIMER IS THE TIMER ROUTINE 
C THAT RETURNS THE CURRENT PROCESSOR TIME ON A SUN/ 3 
C WORKSTATION RUNNING THE UNIX OPERATING SYSTEM AT THE 

C UNIVERSITY OF TENNESSEE KNOXVILLE. 

0 *************************** ********************************************** 

C 

C INPUT PARAMETER - 

C IDUMMY . A DUMMY INTEGER VARIABLE. 

C 

C PROGRAM SUBROUTINE . 

C GTIMBR. 

C 


0 ******************************** ********* ****************** ************* ****** 
c 

REAL FUNCTION CTIMB ( IDUMMY ) 

********************** ***** ***** **** ********* *********************** ********* * 


INTEGER IDUMMY, IPRNTE, IPRNTS, MAXINT 
REAL RATIOL, RATIOS, TIME ,X 

******************************************************************** ********** 


COMMON / CSMSYS / IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 

1 TIME 

***************************************** ********* **** ************** **** ****** 


CALL GTIMER ( X ) 
CTIMB = X - TIME 
TIME = X 
RETURN 
C 

END 
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£******************************** ******* *** ******** ****************************** 
£***************************] *************************** ************************** 
C GETJDP GET NUMBER OF JOINTS ... 

Q**** **************************************************************** ************ 
Q**** **************************************************************** ************ 

c 

C PURPOSE - THIS ROUTINE RETRIEVES THE TOTAL NUMBER OF JOINTS 
C FOR THE PROBLEM TO BE SOLVED. 

C 

C PARAMETERS INITIALIZED - 
C NUMJNT - THE TOTAL NUMBER OF JOINTS. 

C 

C ERROR CODES . 

C 0 - ERROR CODES 

C 1013 - INCORRECT EXECUTION SEQUENCE 

C 1014 - THE NUMBER OF ITEMS AVAILABLE FROM THE RETRIEVED 

C DATASET IS LESS THAN TWO. SEE REMARK. 

C 

C 

C REMARK • 

C THE CURRENT VERSION OF TESTBED DATABASE ASSUMES THAT 

C ALL JOINTS HAVE THE MAXIMUM DEGREES OF FREEDOM, THE 

C NUMBER OF JOINTS AND THE MAXIMUM DEGREES PER JOINT IS 

C FROM THE FIRST TWO ITEMS RETRIEVED. IN CASE OF 

C VARIABLE DEGREES OF FREEDOM PER JOINT, DUMMY DATA IS 

C STORED. 

C 

C 

C PROGRAM SUBROUTINES - 

C QKINFO, GETRECI, BMSG 

C 

C CSM TESTBED DATASETS ACCESSED . 

C JDFl.BTAB.* 

C 

Q** ********************************************************************** ******** 

c 

SUBROUTINE GBTJDF ( IBUF ) 

C 

INTEGER*! IBUF(l) 

C 

c 

Q****************** ****************************************************** ******** 

c 

CHARACTER* 40 LIBNAM 

CHARACTER* 61 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
INTEGER*4 IPRNTE, IPRNTS, MAXINT 

INTEGER*4 MSGLVL, IERR, MAXCSM 

INTEGER*4 BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER*4 MAXDOF, NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 

C 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON / CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /PRBLEM/ MAXDOF, NEQNS .NUMJNT 

******************************************************************************** 
INTEGER* 4 LEN 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 11 ) 

1 FORMAT ( / 5X, ’GETJDF - GET NUMBER OF JOINTS AND ... » ) 

IF (( STAGE .LT. 10 ) .OR. ( IERR .NE. 0 )) GO TO 100 


EACH DATASET IS IDENTIFIED BY A STRING OF 
’MAINKE Y. EX TENSION. CYCLE 1. CYCLE 2. CYCLE 3* 
MAXIMUM NUMBER OF CHARACTERS CONTAINED IS 51 


CALL QKINFO ( JDFSET ) 
IF ( IERR .NE. 0 ) RETURN 

STAGE = 15 


GET THE FIRST TWO ITEMS OF THE FIRST RECORD 
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c 

LBN = 2 

CALL GTRBCI ( 1, IBUF, LBN ) 
IF ( IERR .NE. 0 ) RETURN 
IF ( LBN .LT. 3 ) GO TO 200 
C 

NUMJNT = IBUF(l) 


READ IN MAX UNCONSTRAINED DEGREES OF FREEDOM OF THE MODEL 


MAXDOF = IBUF(2) 
STAGE = 20 
RETURN 


ERROR HANDLING 

100 CONTINUE 

IERR = 1013 

IF ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

200 IERR = 1014 

IF ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

END 
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£ ************* ************** ********* ****************** ****************** ******* 
£1 ******** * **** DC* W** ******* ** ********* **** *************************************** 

C GETDOP GET DEGREES OF FREEDOM 

£ ********* **** ***** **** ***** ***** ***m ********* **** ********** ********* ********* ** 

Q ***« ************** *********************** ********* **** ********* ***** ********* ** 

C 

C PURPOSE - THIS ROUTINE RETRIEVES THE DEGREE OF FREEDOM 
C FOR EACH INDIVIDUAL JOINT FROM THE DATABASE. 

C 

C PARAMETERS INITIALIZED - 

C ID OF - IDOF(K) STORES THE STARTING EQUATION NUMBER FOR 
C JOINT K. THE DEGREES OF FREEDOM FOR JOINT K IS 

C GIVEN BY IDOF(K + l) - IDOF(K). THE TOTAL NUMBER 

C OF EQUATIONS IS EQUAL TO ID OF(NUMJNT + 1 ) - 1. 

C MAXDOF - THE MAXIMUM DEGREE OF FREEDOM RETRIEVED FOR AN 

C INDIVIDUAL JOINT. 

C NEQNS • THE NUMBER OF EQUATIONS EQUALS THE TOTAL DEGREES 

C OF FREEDOM. 

C 

C CSM TESTBED DATASET ACCESSED . 

C CURRENTLY NONE. 

C 

q ************************************************************************* 

C WARNING 

£ ************************************************************************* 

C THIS SUBROUTINE MUST BE MODIFIED FOR PROBLEMS WITH 

C VARIABLE DEGREES OF FREEDOM PER NODE 

C 


q ****************** ************** ********* **** ********* ***** ********* *********** 

C 

SUBROUTINE GETDOF ( IDOF, IBUF ) 

C 

INTEGER*! IDOF(l), IBUF(l) 

CHARACTERS LIBNAM 

CHARACTER*51 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

INTEGERS IPRNTE, IPRNTS, MAXINT 

INTEGER*! MSGLVL, IERR, MAXCSM 

INTEGER** BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER** MAXDOF, NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON / CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /PRBLEM/ MAXDOF, NEQNS , NUMJNT 

********************************************* ********* ************************** 

INTEGER** DEGREE, I 
C 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 11 ) 

11 FORMAT( / 5X, ’GETDOF - GET DEGREES OF FREEDOM ... ' ) 

C 

IF (( STAGE .LT. 20 ) .OR. ( IERR .NE. 0 )) GO TO 500 
C 

c 

C THE FOLLOWING LINES OF CODE IS TEMPORARY 

C FOR THE-FIXED DEGREE PROBLEMS 

C — 

DEGREE = MAXDOF 
IDOF(l) = 1 

DO 100 I = 2, NUMJNT+1 

IDOF(I) = ID OF(I-l) + DEGREE 

IF ( MAXDOF .LT. DEGREE ) MAXDOF = DEGREE 
100 CONTINUE 

NEQNS = IDOF(NUMJNT+l) - 1 
STAGE = 30 
RETURN 
C 

500 CONTINUE 

IERR =1019 

IF ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
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******************************************************************** *********** 
***************************************** **** ************************* 

GTZBRO .... DETECT DUMMY ROWS 

******************************************************************************* 

******************************************************************************* 

PURPOSE . THIS ROUTINE IDENTIFIES DUMMY ROWS (ALL ZEROS) IN 
THE DATA MATRIX. 

INPUT PARAMETERS - 

DOF - AN INTEGER ARRAY OF SIZE EQUAL TO THE TOTAL NUMBER OF 
JOINTS PLUS ONE. 

IDOF(K) STORES THE STARTING EQUATION NUMBER FOR 
JOINT K. THE DEGREES OF FREEDOM FOR JOINT K IS 
GIVEN BY IDOF(K + l) - IDOF(K). THE TOTAL NUMBER 
OF EQUATIONS IS EQUAL TO IDOF(NUMJNT + l) - 1. 

OUTPUT PARAMETERS - 

MASK . THE LINEAR ARRAY MASK STORES A 0 FOR EACH 
ZERO DIAGONAL ELEMENT ENCOUNTERED AND A -1 
FOR EACH NONZERO DIAGONAL ELEMENT. 

WORKING PARAMETERS - 

FBUF . A BUFFER OF MAXIMUM RECORD SIZE FOR RETRIEVING 
REAL OR DOUBLE PRECISION DATA FORM THE TESTBED. 

ERROR CODES . 

1021 - INCORRECT EXECUTION SEQUENCE. 

SUBPROGRAM MODULES - 
QKINFO, GTRECF, EMSG 

CSM TESTBED DATASETS ACCESSED - 
K.SPAR.* 


REMARK - THIS ROUTINE IS NEEDED FOR THE CURRENT RELEASE OF 
TESTBED DATABASE BECAUSE THE CONSTRAINT DATASET DOES NOT 
INCLUDE ZERO ROWS. IN ADDITION, NOTE THAT CURRENTLY 
THE TESTBED STORES MAXDOF EQUATIONS PER JOINT. THEREFORE, 
DUMMY ROWS MUST BE INSERTED FOR THE JOINTS WITH DEGREES 
LESS THAN MAXDOF. 


******************************************************************************* 


SUBROUTINE GTZERO ( DOF, FBUF, MASK ) 

DOUBLE PRECISION FBUF(l) 

INTEGERS MASK(l), DOF(l) 

******************************************************************************* 


CHARACTER* 40 
CHAR A CTER*51 
CHARACTER* 4 


LIBNAM 

JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
RTYPB 


INTEGER*4 IPRNTE, IPRNTS, MAXINT 

INTEGER*4 IDSN , LDI , NLEN , NREC .TRACE 
INTEGER*4 BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER*4 MSGLVL, IERR, MAXCSM 

INTEGER*4 MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 


COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM/ MAXDOF , NEQNS , NUMJNT 


********************************************************************** 


INTBGER*4 CONRNG, I, II, IROW, IS, ITEMS, JGRPS, JOINT, LEN 
INTEGERS CJNT, NROWS, NCOLS, ISIZE, KOUNT, OVERHD, NZEROS 
DOUBLE PRECISION COEF 
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11 

c 

c 

c 

c 

c 


100 

c 

c 

c 

c 

c 


c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


400 


300 

200 

c 

c 

c 

1 

22 

1 

1 


IP ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 11 ) 

FORMAT( /5X, ’GTZERO - DETECT DUMMY ROWS ... ’ ) 

IP ( ( STAGE .LT. 30 ) .OR. ( IERR .NE. 0 ) ) GO TO 500 


INITIALIZE MASK TO BE .1 


DO 100 1=1, NEQNS 
MASK(I) = *1 
CONTINUE 


EACH DATASET IS IDENTIFIED BY A STRING OP 
’MAINKEY.EXTENSION.CYCLE1.CYCLE2. CYCLES’ 
MAXIMUM NUMBER OF CHARACTERS CONTAINED IS 51 


CALL QKINFO ( KSPAR ) 

IP ( IERR .NE. 0 ) RETURN 

OVERHD = 0 
KOUNT = 0 
NZEROS = 0 
TRACE = TRACE + 10 
DO 200 I = 1, NREC 
LEN = NLEN 

CALL GTRECP ( I, FBUP, LBN ) 
IP ( IERR .NE. 0 ) RETURN 


DETERMINE NUMBER OF JOINT GROUPS IN CURRENT RECORD 


JGRPS = PBUP(l) 

ITEMS = 1 

OVERHD = OVERHD + 1 
DO 300 II = 1, JGRPS 

CONRNG = FBUF(ITEMS-H) 

JOINT = PBUP(ITEMS-f 2) 

NROWS = DOF( JOINT+l) - DOF(JOINT) 


COMPUTE THE SIZE OF DATA ITEMS. IN TOTAL 

CONRNG SUBMATRICBS INCLUDING DIAGONAL SUBMATICES 


ISIZE = 0 

DO 350 IS = 1, CONRNG 

CJNT = FBUF(ITBMS + 1+IS) 

NCOLS = DOF(CJNT + l) - DOF(CJNT) 
ISIZE = ISIZE + NROWS*NCOLS 
CONTINUE 

ITEMS = ITEMS + 1 + CONRNG 
OVERHD = OVERHD + 1 + CONRNG 


ACCESS THE DIAGONAL ELEMENTS ON THE DIAGONAL MATRIX 


IROW DOF( JOINT) . 1 
NCOLS = NROWS 
DO 400 IS = 1, NCOLS 

COEF = FBUF(ITEMS+(IS-l)*NROWS+IS) 


A DUMMY ROW IS DETECTED 


IP ( COBP .EQ. 0.0D0 ) THEN 
MASK ( IROW + IS ) := 0 
KOUNT = KOUNT + 1 
ENDIP 
CONTINUE 

ITEMS = ITEMS + ISIZE 
NZEROS = NZEROS + ISIZE 
CONTINUE 
CONTINUE 
STAGE = 40 


PRINT DEBUGGIN DATA ... 


IP ( MSGLVL .GE. 3 ) WRITE ( IPRNTS, 22 ) KOUNT, 
OVERHD, NZEROS 

FORMAT ( 15X, ’NUMBER OF DUMMY ROWS: ’ , 18 

/15X, ’K.SPAR.* INDEX OVERHEAD:’, 18 
/15X, ’K.SPAR.* NONZEROS : ’, 18 ) 
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RETURN 


600 CONTINUE 

ERROR HANDLING ... 

IERR s 1021 

IF ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

END 
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£<**** ********* mm*** ********* ***** ************* ********* ***** **** ***** **** ***** *** 
q************* **************************************************************** *** 

GTCOND .... RETRIEVE CONSTRAINT INFO 
******************************************************************************** 
****************** **** ******************************** ***** ********************* 

PURPOSE . THIS ROUTINE RETRIEVES THE CONSTRAINED COMPONENTS 
OF EACH JOINT AND TREATS THE DUMMY ROWS AS CONSTRAINED 
TO BE ZERO. 

INPUT PARAMETERS . 

DOF . AN INTEGER ARRAY OF SIZE EQUAL TO THE TOTAL NUMBER OF 
JOINTS PLUS ONE. 

IDOF(K) STORES THE STARTING EQUATION NUMBER FOR 
JOINT K. THE DEGREES OF FREEDOM FOR JOINT K IS 
GIVEN BY IDOF(K-f-l) - IDOF(K). THE TOTAL NUMBER 
OF EQUATIONS IS EQUAL TO IDOF(NUMJNT + l ) - 1. 

MASK - RECORD OF DUMMY ROWS. 

OUTPUT PARAMETERS - 

MASK - RECORD OF CONSTRAINED VARIABLES IN ADDITION TO 
DUMMY ONES. 

CSIZE - TOTAL NUMBER OF NONZERO CONSTRAINTS. 

WORKING PARAMETERS - 

1BUF - A BUFFER OF MAXIMUM RECORD SIZE FOR RETRIEVING 
INTEGER DATA FORM THE TESTBED. 

KC . AN TEMPORARY INTEGER ARRAY OF SIZE (MAXDOF+1) 

NEEDED IN DECODING THE CONSTRAINT DATA. 

ERROR CODES - 

1022 - INCORRECT EXECUTION SEQUENCE. 

SUBPROGRAM MODULES . 

QKINFO, GTRECI, DECODE, EMSG 

CSM TESTBED DATASETS ACCESSED - 

CON..* OR CON..i (IF MULTIPLES EXISTS IN DATA LIBRARY) 


REMARKS . 

IT IS ASSUMED THAT THE CONSTRAINED DATA IS STORED 
IN THE DATASET IN THE ORDER OF JOINT NUMBERS. 


************* 


***************************************************************** 


SUBROUTINE GTCOND ( DOF, IBUF, KC, MASK, CSIZE ) 

INTEGERS DOF(l), IBUF(l), KC(l), MASK(l), CSIZE 
******************************************************************************** 


CHARACTER* 40 
CHARACTER* 51 
CHARACTERS 


LIBNAM 

JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
RTYPE 


INTEGER*4 IPRNTE, IPRNTS, MAXINT 
INTEGERS IDSN , LDI , NLEN , NREC , TRACE 

INTEGERS BUFMAX, MXUSBD, MXREQD, STAGE 

INTEGER*4 MSGLVL, IERR, MAX CSM 

INTEGER*4 MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 


COMMON / CSMS YS / IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM/ MAXDOF , NEQNS , NUMJNT 


a******************************)************** ************** ********************* 


INTEGERM I, II, IROW, JOINT, K, LEN, DEGREE, ZKOUNT,FKOUNT, 
1 ZDUMMY 

C 
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IP ( MSGLVL .GB. 2 ) WRITE ( IPRNTS, II ) 

11 FORMAT( / 5X, '6TCOND . GET CONSTRAINTBD VARIABLES... ’ ) 
C 

IP ( ( STAGE .LT. 40 ) .OR. ( IERR .NE. 0 ) ) GO TO 500 
C 

C 

C EACH DATASET IS IDENTIFIED BY A STRING OP 
C , MAINKEY.EXTENSION.CYCLEl.CYCLE2. CYCLES’ 

C MAXIMUM NUMBER OP CHARACTERS CONTAINED IS 61 

CALL QKINFO ( CON ) 

IP ( IERR .NE. 0 ) RETURN 
TRACE = TRACE + 10 

C KOUNTING NONZERO CONSTRAINTS 

C 

CSIZE = 0 

C 

C KOUNTING ZERO CONSTRAINTS 

C 

ZKOUNT = 0 
JOINT = 1 
DO 100 I = 1, NRBC 
LBN = NLEN 

CALL GTRECI ( I, IBUF t LBN ) 

IP ( IERR .NE. 0 ) RETURN 
DO 200 II = 1, LBN 

IP ( JOINT .GT. NUMJNT ) GO TO 200 

C CONSTRAINTS ARB ENCODED INTO 7 BITS 

C WHICH ARE DECODED INTO AN INTEGER 

C ARRAY KC OP SIZE 7 f 

C 

CALL DECODE ( IBUF(II), KC ) 

DEGREE = DOF( JOINT-f 1) - DOF(JOINT) 

IROW = DOF( JOINT) - 1 
DO 300 K = 1, DEGREE 
IP ( KC(K) .EQ. 1 ) THEN 

C 

C ZERO CONSTRAINTS 

MASK(IROW-fK) = 0 
ZKOUNT = ZKOUNT + 1 
ELSE IP ( KC(K) .BQ. 2 ) THEN 

C NONZERO CONSTRAINTS 

MASK(IROW+K) = 1 
CSIZE = CSIZE + 1 
ENDIP 

300 CONTINUE 

JOINT = JOINT + 1 
200 CONTINUE 

100 CONTINUE 

C KOUNTING UNCONSTRAINED DEGREES OF FREEDOM AND 

C THE NET ZERO CONSTRAINTS INCLUDING DUMMY ROWS 
C 

FKOUNT = 0 
ZDUMMY = 0 
DO 400 1 = 1, NEQNS 

IP ( MASK(I) .EQ. -1 ) PKOUNT = FKOUNT + 1 
IF ( MASK(I) .EQ. 0 ) ZDUMMY = ZDUMMY + 1 
400 CONTINUE 
STAGE = 50 

C PRINT DEBUGGING DATA ... 

IF ( MSGLVL .GB. 3 ) WRITE (IPRNTS, 22) ZKOUNT, CSIZE, 

1 FKOUNT, ZDUMMY 

22 FORMAT( 15X, 26H ZERO CONSTRAINTS ARB , 18 
1 /15X, 26HNONZERO CONSTRAINTS ARE , 18 

1 / 1 5X, 26HFRBE VARIABLES ARE , 18 

1 / 15X, 26HDUMMY ROWS + 0 CONSTRAINTS, 18 ) 

RETURN 

C 

500 CONTINUE 


100 


I 


o o o 


ERROR HANDLING 
IERR = 1022 

IP ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

END 
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£ *********** ** ***** **** ***** ***** **** * **** **** ***** **** ***** **** ***** ********* *** 
£******************************************************************************** 
C GTMOTI ....GET NONZERO CONSTRAINTS 

£***************************************** ********* ****************** ************ 
£********************************************************************♦*********** 
c 

C PURPOSE - TO RETRIEVE NUMERIC FOR NONZERO CONSTRAINTS. 

C 

C INPUT PARAMETERS 

C MASK - CONSTRAINT INFORMATION FOR EACH VARIABLE. 

C 

C OUTPUT PARAMETERS 

C MASK • THE LOCATIONS CORRESPONDING TO NONZEROR CONSTRAINTS 

C CONTAIN A POINTER TO THE NUMERIC VALUE IN FCON. 

C FCON - AN ARRAY OF CSIZB FLOATING-POINT CONSTRAINTS. 

C 

c 

C WORKING PARAMETERS 

C FBUF • A REAL OR DOUBLE PRECISION BUFFER OF SIZE BUFMAX. 

C THE ACTUAL TYPE IS AS DECLARED. 

C 

C ERROR CODES - 

C 1025 - INCORRECT EXECUTION SEQUENCE. 

C 1026 - UNEXPECTED NONZERO CONSTRAINT VALUE. 

C 1027 - ZERO ENTRY FOR A NONZERO CONSTRAINT OCCURS. 

C 

C SUBROUTINE PROGRAMS - 

C QKINFO, GTRECF, EMSG. 

C 

C CSM TESTBED DATASETS ACCESSES • 

C APPL.MO Tl.i.j. 

C 

C 

C REMARKS - 

C IT IS ASSUMED THAT THE CONSTRAINT VALUES ARE STORED 

C IN SEQUENCE FROM 1 TO NEQNS. 

C 

c 

£************* *********************** ********* ********* ************** ********* *** 

c 

SUBROUTINE GTMOTI ( FBUF, MASK, FCON, CSIZE ) 

C 

INTEGER** MASK(l), CSIZE 
DOUBLE PRECISION FBUF(l), PCON(l) 

******************************************************************************** 
CHAR A CTERMO LIBNAM 

CHARACTER 41 51 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

CHARACTER* 4 RTYPE 

INTBGERM IPRNTB, IPRNTS, MAXINT 

INTBGERM IDSN , LDI , NLEN , NREC , TRACE 

INTBGERM BUFMAX, MXUSED, MXREQD, STAGE 

INTBGERM MSGLVL, IERR, MAXCSM 

INTBGERM MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM/ MAXDOF , NEQNS , NUMJNT 

****************** ********* ************** ********* ********* ********************* 

INTBGERM NITEMS, KPTR, LEN, I, J 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 11 ) 

11 FORMAT( / 6X, 'GTMOTI . GET NONZERO CONSTRAINTS... ’ ) 

IF (( STAGE .LT. 50 ) .OR. ( IERR .NE. 0 ) ) GO TO 1000 

IF ( CSIZE .EQ. 0 ) THEN 


NONZERO CONSTRAINTS ARE NOT EXPECTED 


102 



c 


IP ( MSGLVL .GE. 3 ) WRITE ( IPRNTS, 21 ) 

21 FORMAT( /10X, ’APPLIED DISPLACEMENTS ARE NOT EXPECTED.’) 

STAGE = 60 
RETURN 
ENDIP 
C 

C RETRIEVE NEQNS ITEMS FORM ’APPL.MOTI * * 

CALL QKINFO ( APPLM ) 

IP ( IERR .NE. 0 ) RETURN 

TRACE = TRACE + 10 

NITEMS = 0 

KPTR = 0 

DO 100 I = 1, NREC 

LEN = MINO ( NEQNS - NITEMS, NLEN ) 

IP ( LEN .GT. 0 ) THEN 

CALL GTRECP ( I, FBUF, LEN ) 

IP ( IERR .NE. 0 ) RETURN 
DO 200 J = 1, LEN 

NITEMS = NITEMS + 1 

C 

C CHECK ERROR DUE TO INCONSISTENT CONSTRAINT VALUES 

C 

IP (( MASK(NITEMS) -NE. 1 ) .AND. 

1 ( FBUF(J) .NE. 0.0D0 )) GO TO 1100 

IP (( MASK (NITEMS) .EQ. 1 ) .AND. 

1 ( FBUF(J) .EQ. 0.0D0 )) GO TO 1200 

IP ( MASK (NITEMS) .EQ. 1 ) THEN 
C 

C ENTER NUMERIC FOR NONZERO CONSTRAINT 

KPTR = KPTR + 1 

FCON(KPTR) ss PBUP(J) 

C STORE THE ADDRESS POINTER IN MASK 

MASK(NITEMS) = KPTR 
ENDIP 

200 CONTINUE 

ENDIP 

100 CONTINUE 
STAGE s 60 
RETURN 
C 

C ERROR HANDLING 

1000 CONTINUE 
IERR = 1025 

IP ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

1100 CONTINUE 
IERR = 1026 

IP ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

1200 CONTINUE 
IERR = 1027 

IP ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

END 
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********* ********* **** *** * * **** * * * ** * ** **** ** ***** **** ** *** ***** * ************* * 
******************************************************************************* 

GETIJ .... INPUT NONZERO STRUCTURES 
******************************************************************************* 
******************************************************************************* 


PURPOSE . TO RETRIEVE NONZERO STRUCTURES FROM DATASET KMAP..* 
AND INPUT THE SAME TO SPARSPAK-A. 

INPUT PARAMETERS 

DOF - AN INTEGER ARRAY OF SIZE EQUAL TO THE TOTAL NUMBER OF 
JOINTS PLUS ONE. 

IDOF(K) STORES THE STARTING EQUATION NUMBER FOR 
JOINT K. THE DEGREES OF FREEDOM FOR JOINT K IS 
GIVEN BY IDOF(K-j-l) - IDOF(K). THE TOTAL NUMBER 
OF EQUATIONS IS EQUAL TO IDOF(NUMJNT-J-l) - 1. 

MASK - CONSTRAINT INFORMATION FOR EACH VARIABLE. 

OUTPUT PARAMETERS 

S - NONZERO STRUCTURES SET UP BY SPARSPAK-A. 

WORKING PARAMETERS 

IBUF - AN INTEGER BUFFER OF SIZE BUFMAX. 

ICLQ - A TEMPORARY ARRAY OF SIZE MAXDOF. 

ERROR CODES - 

1023 - INCORRECT EXECUTION SEQUENCE. 

SUBROUTINE PROGRAMS - 
QKINFO, GTRBCI, EMSG 

SPRSPAK-A SUBROUTINES . 

IJBEGN, INCLQ, INIJ, IJEND. 

CSM TESTBED DATASETS ACCESSES . 

KMAP. * 


K* **************************************************************** ************ 

SUBROUTINE GETIJ ( DOF, IBUF, ICLQ, MASK, S ) 

INTEGER** DOF(l), IBUF(l), ICLQ(l), MASK(l), S(l) 
****************************************************************************** 


CHARACTER* 40 
CHARACTER* 51 
CHARACTER* 4 


LIBNAM 

JDFSBT, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
RTYPE 


INTBGER*4 IPRNTE, IPRNTS, MAXINT 
INTBGER*4 IDSN , LDI , NLEN , NREC , TRACE 

INTEGER*4 BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER*4 MSGLVL, IERR, MAXCSM 

INTEGER*4 MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 


COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK / IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSBT, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM/ MAXDOF , NEQNS , NUMJNT 


:********************************************************** ****************** *** 


INTEGER*4 CONRNG, I, II, ICOL, IROW, ITEMS, J, JGRPS, JOINT, 
1 K, NCLQ, LBN, IX, JX, LRNG, NODES, JJ, NROWS, 

1 NCOLS 

IF ( MSGLVL .GB. 2 ) WRITE ( IPRNTS, 11 ) 

FORMAT( / 5X, ’GETIJ - INPUT NONZERO STRUCTURES... ’ ) 

IF (( STAGE .LT. 60 ) .OR. ( IERR .NB. 0 ) ) GO TO 1000 

CALL IJBEGN 


INIJ INSURES NONZERO FOR ALL DIAGONAL ELEMENTS 
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IF POSITION (NEQNS, NEQNS) IS ENTERED 


CALL INIJ ( NEQNS, NEQNS , S ) 


ACCESS EACH RECORD IN DATA SET 'KMAP..* • 


CALL QKINFO ( KMAP ) 

IF ( IBRR .NE. 0 ) RETURN 
TRACE = TRACE + 10 
DO 100 I = 1, NREC 
LEN = NLEN 

CALL GTRECI ( I, IBUF, LEN ) 
IF ( IERR .NE. 0 ) RETURN 


DETERMINE NUMBER OF JOINT GROUPS IN CURRENT RECORD 

JGRPS = IBUF(l) 

ITEMS = 1 

DO 200 II = 1, JGRPS 


GET THE CURRENT JOINT AND COMPUTE THE ROW NUMBER 


JOINT = IBUF(ITEMS-fl) 

NUMBER OF DEGREES FOR CURRENT JOINT 


NROWS as D OF( JOINT + l) - DOF(JOINT) 

COMPUTE THE THE ROW NUMBER BY IROW + K, 
WHERE IROW IS GREATER THAN OR EQUAL TO 0 


IROW = DOF( JOINT) - 1 
NCLQ = 0 

DO 300 K = 1, NROWS 

IF ( MASK ( IROW + K ) .EQ. -1 ) THEN 


THIS ROW IS NOT CONSTRAINED 


NCLQ = NCLQ + 1 
ICLQ(NCLQ) = IROW + K 
ENDIF 
CONTINUE 


INPUT DIAGONAL BLOCK TO SPARSPAK 


IF ( NCLQ .GT. 0 ) CALL INCLQ( NCLQ, ICLQ, S ) 

SKIP UNRELATED ITEMS IN CURRENT JOINT GROUP 


LRNG = IBUF(ITEMS + 2) 

ITEMS = ITEMS + 2 
DO 220 JJ = 1, LRNG 

NODES = IB UF( ITEMS + 1) 

ITEMS - ITEMS + 6 + (NODES*(NODES + l))/2 
CONTINUE 

NUMBER OF SUBMATRICES FOR THE CURRENT JOINT 

CONRNG = IBUF(ITEMS + 1 ) 

ITEMS = ITEMS + 1 


ENTER NONZERO IN THE CONNECTED SUBMATRIX 
IN ADDITION TO THE DIAGONAL SUBMATRIX 


DO 400 J = 1, CONRNG. 1 
JOINT = IBUF(ITEMS + J) 


DEGREE OF FREEDOM OF THE CONNECTED JOINT 
NCOLS = DOF( JOINT -fl). DOF(JOINT) 

COMPUTE STARTING COLUMN NUMBER 
ICOL = DOF( JOINT) - 1 

COMPUTE NONZERO POSITION COLUMN BY COLUMN 
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DO 500 JX = 1, NCOLS 
DO 550 IX = 1, NROWS 

IF ( ( MASK(ICOL+JX) .EQ. * I ) .AND. 

1 ( MASK(IROW + IX) .EQ. -1 ) ) THEN 

C 

C THE CORRESPONDING VARIABLES 

C ARE NOT CONSTRAINED 

C 

CALL INIJ ( IROW+IX, ICOL-f JX, S ) 
ENDIF 

550 CONTINUE 

600 CONTINUE 

400 CONTINUE 

ITEMS = ITEMS + 2*CONRNG - 1 

C 

C END OF CURRENT JOINT GROUP 

C 

200 CONTINUE 

C 

C END OF CURRENT RECORD 

C 

100 CONTINUE 
C 

CALL IJBND ( S ) 

STAGE = 70 
RETURN 
C 

1000 CONTINUE 
C ERROR HANDLING 

IERR = 1023 

IF ( MSGLVL .GB. 2 ) CALL BMSG 
RETURN 
C 

END 
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:************ ***** ********* ***** ********* ********* ********* ********* ********* *** 
:**«* ** *** Dc next ******** >K ********** ************* ** ******* ***** ** ** ************ * * ** * 

GTFORC .... INPUT RIGHT HAND SIDE 

:******************************************************************************* 
[***************** *********************************************************** *** 

PURPOSE - TO RETIE VE RIGHT HAND SIDE FROM DATASET APPL.FORC.i.j 
AND INPUT THOSE COMPONENTS CORRESPONDING TO UNCONSTRAINED 
VARIABLES TO SPARSPAK-A. 

INPUT PARAMETERS 

MASK - CONSTRAINT INFORMATION FOR EACH VARIABLE. 

S - INPUT TO SPARSPAK-A ROUTINES. 

OUTPUT PARAMETER 

S - SPARSPAK-A OUPUT. 

WORKING PARAMETERS 

FBUF - A REAL OR DOUBLE PRECISION BUFFER OF SIZE BUFMAX. 

THE ACTUAL TYPE IS AS DECLARED. 

ERROR CODES - 

1024 - INCORRECT EXECUTION SEQUENCE. 

SUBROUTINE PROGRAMS . 

QKINFO, GTRECF, EMSG. 

SPRSPAK-A SUBROUTINES - 
INBI. 

CSM TESTBED DATASETS ACCESSES • 

APPL.FORC.i.j. 


REMARKS . 

IT IS ASSUMED THAT THE ROWS CORRESPONDING TO DUMMY AND 
CONSTRAINED VARAIBLES ARE INCLUDED IN THE DATA MATRIX. 


t******** *********************** ****************** ****************************** 

SUBROUTINE GTFORC ( FBUF, MASK, S ) 

INTEGER** MASK(l) 

DOUBLE PRECISION FBUF(l), S(l) 

i********************* ************************************************** ******** 


CHARACTERMO LIBNAM 

CHARACTER’ 1 ' 51 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
CHARACTER** RTYPE 

INTEGER** IPRNTE, IPRNTS, MAXINT 

INTEGER** IDSN , LDI , NLEN , NREC .TRACE 
INTEGER*! BUFMAX, MXUSED, MXRBQD, STAGE 

INTEGER*! MSGLVL, IERR, MAXCSM 

INTEGER** MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 


COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK / IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM/ MAXDOF , NEQNS , NUMJNT 

[************ ******************************************************* ************ 


INTEGER** I, J, IROWS, LEN 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 11 ) 

11 FORMAT( / 5X, 'GTFORC - INPUT RIGHT HAND SIDE... • ) 

IF {( STAGE .LT. TO ) .OR. ( IERR .NE. 0 ) ) GO TO 1000 


RETRIEVE RIGHT HAND SIDE FORM 'APPL.FORC.* * 
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CALL QKINFO ( APPLF ) 

C NOTE APPLY .FORC., DOES NOT NECESSARILY EXIST 

IF ( IERR .NE. 0 ) GO TO 900 
TRACE m TRACE + 10 
IROWS = 0 
DO 100 1=1, NREC 

LEN = MIN0 ( NBQNS - IROWS, NLEN ) 

IF ( LBN ,GT. 0 ) THEN 

C 

C READ NEXT RECORD 

C 

CALL GTRECF ( I, FBUF, LEN ) 

IF ( IERR .NE. 0 ) RETURN 

C RETRIEVE EACH ITEM IN CURRENT RECORD 

C 

DO 200 J = 1, LEN 
IROWS = IROWS + 1 
IF ( MASK ( IROWS ) .EQ. -1 ) THEN 

C 

C THE VARIABLE IS NOT CONSTRAINED 

CALL INBI ( IROWS, FBUF(J), S ) 

ENDIF 

200 CONTINUE 

ENDIF 

100 CONTINUE 
STAGE = 80 
RETURN 
C 

900 CONTINUE 

C 

C RIGHTHAND SIDE DOES NOT EXIST 

IF ( MSGLVL .GE. 3 ) WRITE ( IPRNTS, 21 ) 

21 FORMAT( /1DX, ’THERE IS NO APPLIED FORCE VECTOR’) 
IERR = 0 
STAGE = 80 
RETURN 
C 

1000 CONTINUE 
C ERROR HANDLING 

IERR = 1024 

IF ( MSGLVL .GB. 2 ) CALL EMSG 
RETURN 
C 

END 
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******************************************************************************* 
****************** **** ** *** He**** **** ***** ********* **** ***** * ******** ********* ** 

GTNUM5 ... INPUT NONZERO NUMERICS 
******************************************************************************* 
******************************************************************************* 


PURPOSE - TO RETRIEVE AND INPUT NUMERICAL NONZEROS OP THE 
SYSTEM MATRIX. IN ADDITION, RIGHT HAND SIDE IS APPROP- 
RIATELY ADJUSTED USING CONSTRAINTS AVAILABLE. 

INPUT PARAMETERS 

DOF . AN INTEGER ARRAY OF SIZE EQUAL TO THE TOTAL NUMBER OF 
JOINTS PLUS ONE. 

IDOF(K) STORES THE STARTING EQUATION NUMBER FOR 
JOINT K, THE DEGREES OF FREEDOM FOR JOINT K IS 
GIVEN BY IDOF(K + l) - IDOF(K). THE TOTAL NUMBER 
OF EQUATIONS IS EQUAL TO IDOF(NUMJNT+l) - 1. 

MASK - THE LOCATIONS CORRESPONDING TO NONZEROR CONSTRAINTS 
CONTAIN A POINTER TO THE NUMERIC VALUE IN FCON. 

THE OTHER LOCATIONS INDICATE FREE OR CONSTRAINED 
TO ZERO VARIABLES. 

FCON - AN ARRAY OF CSIZE FLOATING-POINT CONSTRAINTS. 

S - STORAGE ARRAY FOR SPARSPAK-A. 

WORKING PARAMETERS 

FBUF - A REAL OR DOUBLE PRECISION BUFFER OF SIZE BUFMAX. 

THE ACTUAL TYPE IS AS DECLARED. 

ERROR CODES - 

1028 - INCORRECT EXECUTION SEQUENCE. 

SUBROUTINE PROGRAMS - 
QKINFO, GTRBCF, EMSG. 

SPARSPAK-A ROUTINES - 
INAIJ5, INBI. 

CSM TESTBED DATASETS ACCESSES - 
K.SPAR.*. 


REMARKS - 

IT IS ASSUMED THAT THE VARIABLES ARE ORDERED IN THE 
GIVEN ORDER OF THE JOINTS AND DEGREES. 


******************************************************************************** 


SUBROUTINE GTNUM5 ( DOF, FBUF, MASK, FCON, S ) 

INTEGERS DOF(l), MASK(l) 

DOUBLE PRECISION FBUF(l), FCON(l), S(l) 

******************************************************************************** 


CHARACTER* 40 
CHARACTER* 61 
CHARACTER* 4 


LIBNAM 

JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
RTYPB 


INTBGER*4 IPRNTB, IPRNTS, MAXINT 
INTEGER*4 IDSN , LDI , NLEN , NREC , TRACE 

INTEGER*4 BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER*4 MSGLVL, IERR, MAXCSM 

INTEGER*4 MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 


COMMON /CSMSYS/ IPRNTB, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM / MAXDOF , NEQNS , NUMJNT 


******************************************************************************* 


INTEGER*! CONRNG, I, II, ICOL, IROW, ISTRT, ITEMS, 
1 JGRPS, JOINT, M, MTXKNT, MYI, MYJ, NCOL, 
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c 


NROW, LBN , NCOLS, NROWS 
DOUBLE PRECISION COEP, BIX, BJX 


IP ( MSGLVL .GB. 2 ) WRITE ( IPRNTS, 11 ) 

11 FORMAT( / 5X, ’GTNUM5 - GET NONZERO NUMERIC... ' ) 

C 

IP (( STAGE .LT. 80 ) .OR. ( IERR .NE. 0 ) ) GO TO 1000 
C 

C 

C ACCESS EACH RECORD IN DATA SET ’K.SPAR.* * 

C 

CALL QKINFO ( KSPAR ) 

IP ( IERR .NE. 0 ) RETURN 
TRACE = TRACE + 10 
DO 100 I = 1, NREC 
LEN = NLEN 

CALL GTRBCP ( I, PBUP, LEN ) 

IP ( IERR .NE. 0 ) RETURN 


DETERMINE NUMBER OP JOINT GROUPS IN CURRENT RECORD 


JGRPS = FBUF(l) 
ITEMS = 1 

DO 200 II = 1, JGRPS 


GET NUMBER OP SUBMATRICES 


CONRNG = FBUF(ITEMS + 1) 


GET THE CURRENT JOINT 


JOINT = FBUP(ITEMS+2) 

IROW = DOF( JOINT) - 1 

NROWS = DOF(JOINT + l) - DOF(JOINT) 

ISTRT = ITEMS + 1 + CONRNG 


RETRIEVE UPPER TRIANGULAR PART OP DIAGONAL SUBMATRIX 


NCOLS = NROWS 
DO 400 NCOL = 1, NCOLS 
MYJ = IROW + NCOL 
DO 500 NROW = 1, NROWS 
ISTRT = ISTRT + 1 
IP ( NROW .GT. NCOL ) GO TO 500 
COBP = PBUP(ISTRT) 

MYI = IROW + NROW 


RETRIEVE THE NONZERO CONSTRAINTS 


IP ( MASK(MYI) .GT. 0 ) BIX = FCON(MASK(MYI)) 
IP ( MASK(MYJ) .GT. 0 ) BJX = FCON(MASK(MYJ)) 
IP ( MYI .EQ. MYJ ) THEN 

IF ( MASK(MYI) .NE. .1 ) THEN 


CHANGE DIAGONAL ELEMENT TO BE 1.0D0 
FOR CONSTRAINED ROW 


COEP = 1.0D0 


ENTER NONZERO CONSTRAINT VALUE AS RHS 


IP ( MASK(MYI) .GT. 0 ) 

1 CALL INBI ( MYI, BIX, S ) 

ENDIP 


INPUT DIAGONAL ELEMENT COEP 


CALL INAIJ6 ( MYI, MYI, COEP, S ) 
ELSE IP ((MASK(MYJ) .GT. 0 ) .AND. 

1 (MASK(MYI) .EQ. -1 )) THEN 

CALL INBI ( MYI, -COEF*BJX, S ) 
ELSE IF ((MASK(MYI) .GT. 0 ) .AND. 

1 (MASK(MYJ) .EQ. -1 )) THEN 

CALL INBI ( MYJ, -COEF*BIX, S ) 
ELSE IP ((MASK(MYI) .EQ. -1) .AND. 

1 (MASK(MYJ) .EQ. -1)) THEN 

C 
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C INPUT COBF IN LOWER TRIANGULAR MATRIX 

CALL INAIJ5 ( MYJ, MYI, COEF, S ) 

ENDIF 

600 CONTINUE 

C 

C NEXT COLUMN IN DIAGONAL SUBMATRIX 

400 CONTINUE 

C RETRIEVE OFF-DIAGONAL SUBMATRICES IN THE UPPER 

C TRIANGULAR PART OF THE SYSTEM STIFFNESS MATRIX 

C 

MTXKNT = CONRNG - 1 
IF ( MTXKNT .EQ. 0 ) GO TO 199 
ITEMS = ITEMS + 2 
DO 600 M = 1, MTXKNT 

JOINT = FBUF(ITEMS + M) 

ICOL = DOF( JOINT) . 1 

NCOLS = DOF( JOINT+l) - DOF(JOINT) 

DO 800 NCOL = 1, NCOLS 
MYJ = ICOL + NCOL 
DO 900 NROW = 1, NROWS 
ISTRT = ISTRT + 1 
COEF = FBUF(ISTRT) 

MYI = IROW + NROW 

C RETRIEVE NONZERO CONSTRAINTS 

C 

IF ( MASK(MYI) .GT. 0) BIX = FCON(MA SK(MYI)) 

IF ( MASK(MYJ) .GT. 0) BJX = FCON(MA SK(MYJ)) 

C INPUT COEF OR MODIFY RIGHT HAND SIDE 

IF (( MASK(MYI) .E Q. -1 ) .AND. 

1 ( MASK(MYJ) .ECJ. -1 )) THEN 

C 

C ENTER COEF WITH SYMMETRIC POSITION 

C IN LOWER TRIANGULAR TO SPARSPAK-A 

C 

IF ( MYI .LT. MYJ ) 

1 CALL INAIJ5 ( MYJ, MYI, COEF, S ) 

IF ( MYI .GT. MYJ ) 

1 CALL INAIJ5 ( MYI, MYJ, COEF, S ) 

ELSE IF ( (MASK(MYI) .GT. 0 ) .AND. 

1 (MASK(MYJ) .EQ. -1) ) THEN 

CALL INBI ( MYJ, -COEF*BIX, S ) 

ELSE IF ( (MASK(MYJ) .GT. 0 ) .AND. 

1 (MASK(MYI) .EQ. -1) ) THEN 

CALL INBI ( MYI, -COEF*BJX, S ) 

ENDIF 

900 CONTINUE 

C NEXT COLUMN 

800 CONTINUE 

C NEXT SUBMATRIX 

600 CONTINUE 

C PROCESS THE NEXT JOINT GROUP IN THE CURRENT RECORD 

C 

199 ITEMS = ISTRT 

200 CONTINUE 

C NEXT RECORD 

100 CONTINUE 
STAGE = 90 
RETURN 
C 

1000 CONTINUE 
C ERROR HANDLING 
IERR = 1028 
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IP ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

END 
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0 ********************************* *********************************************** 
0 **** **************************************************************** ************ 
C SPACE .... CHECK AVAILABLE STORAGE 

0 ***» ************************************* ********* ****************************** 
0 ************* ******************************************************************* 
c 

C PURPOSE . CHECK STORAGE REQUIRED AGAINST STORAGE AVAILABLE. 

C 

C SUBROUTINE PROGRAMS . 

C EMSG. 

C 

0 ********* ********* ************************************************** **** ******** 
C 

INTEGER FUNCTION SPACE ( IDUMMY ) 

C 

INTEGERS IDUMMY 
C 

CHARACTERS LIBNAM 
CHARACTER'S 1 CDUMMY 
INTEGER** MSGLVL , IERR , MAXCSM 

INTEGER** BUFMAX, MXUSED, MXREQD, STAGE 

C 

COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 
COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
C 

0 ******************************** *************************** ********************* 

c 

IF ( MXREQD .LE. MAXCSM ) THEN 
SPACE = 0 
RETURN 
ELSE 

SPACE = 1 
GO TO 100 
ENDIF 
C 

C ERROR HANDLING 

C 

100 CONTINUE 
IERR = 1001 

IF ( MSGLVL .GE. 2 ) CALL EMSG 
RETURN 
C 

END 
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t ** *** ****** ****** ************** **** ***** ******************** * ***************** * 

LIBOPN .... OPEN DATA LIBRARY ... 

[*** ***** **** ********* **** 1(1 ************** ********* *************************** *** 
iw«* ********* ******************************************************************* 

PURPOSE - THIS ROUTINE OPENS AN EXISTING LIBRARY RESIDENT 
ON A DISKFILB OR MAIN STORAGE, AND CONNECTS IT TO A 
LOGICAL DEVICE INDEX (LDI). THE NAME OP THE LIBRARY 
IS SPECIFIED BY PARAMETER LIBNAM, 

PARAMETERS INITIALIZED - 

LDI - LOGICAL DEVICE INDEX ASSIGNED TO THE EXTERNAL 
DEVICE SPECIFIED BY LIBNAM. 

ERROR CODES - 
0 - NO ERROR. 

1011 - UNSUCCESSFUL OPEN. 

1012 - THE LOGICAL DEVICE NUMBER EXCEEDS THE MAXIMUM VALUE 
OF 30. 

GAL-PROCESSOR ENTRY POINTS . 

LMOPBN, BMSG. 

,*** ******* ****************************** **** 4 < ** * ************* ***** ************ * 


SUBROUTINE LIBOPN 

:************ ***** «»«*** ***** *** *** ****** **** *** * * * * ** ***** ********* **** ******* * 


CHARACTBRMO LIBNAM 

CHAR A CTBR*51 JDFSBT, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

CHARACTER* 4 RTYPE 

INTEGERS IPRNTE, IPRNTS, MAXINT 

INTEGER*4 MSGLVL, IERR , MAXCSM 

INTEGERS IDSN , LDI , NLEN , NREC , TRACE 

INTEGER*4 ICPAD , STAGE 

REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, 
l JDFSBT, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ ICPAD(3), STAGE 

INTBGER*4 LMOPBN 

******* ************************************** **** mm* **** *** ****************** *** 


CHARACTER* 10 LIBKEY 
INTEGER*4 LIMIT 

IF ( MSGLVL .GB. 2 ) WRITE ( IPRNTS, 11 ) LIBNAM 
11 FORMAT ( / 6X, ’LIBOPN- OPEN ’, A40 ) 

IERR = 0 


LIBKEY IS A STRING OF FORM ’MAINKBY/ QUALIFIER* 
MAXIMUM NUMBER OF CHARACTERS IS 10 


LIBKEY = ’ROLD » 

LIMIT = 0 
TRACE = 1000 

LDI = LMOPBN ( LIBKEY, 0, LIBNAM, LIMIT, TRACE ) 


LDI RANGES FROM 1 THROUGH 30 FOR SUCCESSUL OPEN 


IF (( LDI .LT. 1 ) .OR. ( LDI .GT. 30 )) GO TO 100 
STAGE =10 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 21 ) JDFSBT, KMAP, 
1 KSPAR, CON, APPLF, APPLM, STATD 

21 FORMAT(/5X, 35HDATASETS TO BE ACCESSED: 

1 / 5X, 35H 

1 /10X, A51, 

1 /10X, AM, 

1 / 10X, A51, 
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1 /10X, A51, 

1 /10X, AS 1, 

1 /10X, A51, 

1 /10X, AS 1 ) 

RETURN 

C 

100 CONTINUE 

C ERROR HANDLING 

IP ( LDI .LE. 0 ) IERR = 1011 
IP ( LDI .GT. 30) IERR = 1012 
IP ( MSGLVL ,GE, 2 ) CALL EMSG 
RETURN 
C 

END 
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0 ************* ******************* **** ****** ********* ***************************** 
c ******************************************************************** ************ 

C QKINFO ... ANQUIRE DATASET ATTRIBUTES 

0*************************** ************************************************** *** 
0**** ********* ********* ******************************************************* *** 

C 

C PURPOSE - ACQUIRE THE ATTRIBUTES OF A NAMED DATA SET. 

C 

C INPUT PARAMETER - 
C DSNAMB - NAME OF THE DATASET. 

C 

C PARAMETERS UPDATED * 

C IDNS - UNIQUE SEQUENCE NUMBER OF NAMED DATASET. 

C NLEN . LOGICAL LENGTH (ITEMS) OF A RECORD. 

C RTYPE . RECORD TYPE. 

C NREC • TOTAL NUMBER OF RECORDS IN THE DATASET. 

C 

C ERROR CODES . 

C 0 . NO ERROR. 

C 2001 - DATASET DOES NOT EXIST. 

C 2002 - NO RECORD EXISTS IN DATASET. 

C 2003 - RECORD GROUP KEY IS UNDEFINED. 

C 2004 - SEGMENTED RECORD GROUP NOTED. 

C 2000 - RECORD LENGTH GREATER THAN BUFFER LENGTH 
C 

C GAL-PROCESSOR ENTRY POINTS - 
C LMFIND, GMGEKA, GMGECY, BMSG. 

C 

Q****************** ************** **** ******************************************** 

C 

SUBROUTINE QKINFO ( DSNAME ) 

C 

CHARACTER* 51 DSNAME 
C 

0 ***# ******************************************************* ********************* 

c 

CHARACTER*40 LIBNAM 
CHARACTER* 51 CDUMMY 
CHARACTER* 4 RTYPE 
INTEGER*4 MSGLVL, IERR, MAXCSM 
INTEGER*4 IDSN , LDI , NLEN , NREC , TRACE 
INTBGER*4 BUFMAX, MXUSED, MXREQD, STAGE 
C 

COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR, MAXCSM, CDUMMY(7) 
COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
C 

INTEGER*4 LMFIND 
C 

Q************* **************************************************************** *** 

c 

CHARACTER* 1 OP 
CHARACTER* 12 RKEY 
C 

INTEGER* 4 IHI , ILO , MD1M 

C 

C 

C OBTAIN THE SEQUENCE NUMBER OF DATASET DSNAMB 

C MAXIMUM LENGTH OF DSNAMB IS 51 CHARACTERS 

C 

TRACE = TRACE + 10 

IDSN = LMFIND ( LDI, DSNAME, TRACE ) 

IF ( IDSN .EQ. 0 ) GO TO 100 
C 

C OP IS PRESENTLY A DUMMY ARGUMENT FOR BOTH 

C GMGEKA AND GMGECY. 

C 

OP = * ’ 

C 

C RKEY CONTAINS THE RECORD KEY LIFT JUSTIFIED. 

C MAXIMUM LENGTH IS 12 CHARACTERS. 

C 

RKEY s ’DATA ' 

TRACE = TRACE + 10 

C 
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100 


c 

200 


c 

300 


c 

400 


c 

500 


c 


RETRIEVE ATTRIBUTES RTYPE AND NLEN FOR RECORDS OF GIVEN KEY 


CALL GMGEKA ( OP, LDI, IDSN, RKE Y, RTYPE, NLEN, MDIM, TRACE ) 

IF ( NLEN .BQ. 0 ) GO TO 200 

IF ( NLEN .GT. BUFMAX ) GO TO 500 


NUMBER OF RECORDS FOUND WITH GIVEN KEY 


TRACE = TRACE + 10 

CALL GMGECY ( OP, LDI, IDSN, RKBY, NREC, ILO, IHI, TRACE ) 
IF ( NREC .EQ. 0 ) GO TO 300 


NREC = IHI-ILO + 1 FOR AN UNSEGMENTED RECORD GROUP 


IF ( NREC .NE. (IHMLO + 1) ) GO TO 400 
RETURN 


ERROR HANDLING 


CONTINUE 
IERR = 2001 

IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 

CONTINUE 
IERR = 2002 

IF ( MSGLVL .GB. 3 ) CALL EMSG 
RETURN 

CONTINUE 
IERR = 2003 

IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 

CONTINUE 
IERR = 2004 

IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 

CONTINUE 
IERR - 2009 
BUFMAX = NLEN 
IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 

END 
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lew*** **************************************************: 


******************* ************* **** ********* ********* ***** ******************** 

GTRECI ... READ A RECORD PROM A DATASET 
********************** ************** *«-*** **** ********************************** 
****************** ********* *************************** ***** ******************** 


PURPOSE - THIS ROUTINE READS A RECORD FROM A NAMED DATASET. 
THE DATASET MUST BE OF TYPE INTEGER. 


INPUT PARAMETERS - 

RECNUM - RECORD CYCLE OF AN INDIVIDUAL RECORD. 

OUTPUT PARAMETERS. 

LEN • THE NUMBER OF ITEMS CONTAINED IN THE RECORD. 
WORKING PARAMETERS - 

IBUF - A BUFFER OF MAXIMUM RECORD SIZE FOR READIN DATASETS 
OF TYPE INTEGER. 


ERROR CODES - 
0 - NO ERROR. 

2005 - RECORD TYPE IN THE DATASET IS NOT INTEGER. 

2006 - ERROR IN GMGETN DETECTED BY LMERCD. 


GAL-PROCESSOR ENTRY POINTS - 
GMCORN, GMGETN, LMERCD, EMSG. 

******************************************************************************* 


SUBROUTINE GTRECI ( RECNUM, IBUF, LEN ) 

INTEGERS RECNUM, IBUF(l), LEN 

*********** **************************** ********* ********* ******************** 


CHARACTER* 40 DUMMY1 
CHARACTER* 61 CDUMMY 
CHARACTER* 4 RTYPE 

INTEGER*4 IDSN , LDI , NLEN , NREC , TRACE 
INTEGER*4 MSGLVL, IERR, DUMMY2 

COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON /CSMUSR/ DUMMY1, MSGLVL, IERR, DUMMY2, CDUMMY(T) 
INTEGER*4 LMERCD 

:*********** ******************* *********************************************** 


CHARACTERS BUFTYP 
CHARACTER* 12 OP, RKEY 
CHARACTER*24 RNAME 

INTEGER*4 IERROR, IGAP , IHI , ILO , IOFF , MDIM 


DETECT TYPE MISMATCH 


IF ( RTYPE .NE. »I » ) GO TO 500 


CONSTRUCT NAME *RKE Y. RECNUM :RECNUM* FOR AN INDIVIDUAL RECORD 
MAXIMUM LENGTH IS 24 CHARACTERS: 12 FOR RKEY, 5 FOR EACH 
RECNUM REPRESENTING HIGH AND LOW CYCLES. 


RKEY = ’DATA ’ 

ILO = RECNUM 
IHI = RECNUM 

CALL GMCORN ( RNAME, RKEY, ILO, IHI ) 


OP ARGUMENT FOR GMGETx: ’MAINKEY/ QUALIFIER 1 
MAXIMUM LENGTH IS 11: 4 FOR KEY AND 6 FOR QUALIFIER 


OP = ’READ/LENGTH • 

BUFTYP s ’I * 

IGAP = 0 
IOFF = 0 

CALL GMGETN ( OP, LDI, IDSN, RNAME, BUFTYP, IBUF, LEN, MDIM, 
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1 IGAP, IOFF, TRACE ) 


TEST ERROR CONDITION AFTER AN ERROR. SENSITIVE REFERENCE 
TO THE I/O MANAGER 


IERROR = LMERCD ( IERROR ) 
IF ( IERROR ,NB. 0 ) GO TO 600 
RETURN 
C 

500 CONTINUE 
IERR = 2005 

IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 
C 

600 CONTINUE 
IERR = 2006 

IF ( MSGLVL .GB. 3 ) CALL EMSG 
RETURN 
C 

END 
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£************* ***** **** ***** ***** **** ********* ************** ********** * **** ***** 
£************* ***** **** ******************* ************* ************************* 
C GTRECF ... READ A RECORD OP TYPE REAL ... 

£***************************************** ************* ************************* 
£*****«**************** ********** I**!************************#******************** 

c 

C PURPOSE - THIS ROUTINE READS A RECORD PROM A NAMED DATASET. 

C THE DATASET MUST BE OF TYPE REAL OR DOUBLE PRECISION. 

C 

C INPUT PARAMETERS - 

C RECNUM - RECORD CYCLE OF AN INDIVIDUAL RECORD. 

C 

C OUTPUT PARAMETERS. 

C LEN - THE NUMBER OF ITEMS CONTAINED IN THE RECORD. 

C 

C WORKING PARAMETERS . 

C FBUF - A BUFFER OF MAXIMUM RECORD SIZE FOR READIN DATASETS 
C OF TYPE REAL OR DOUBLE PRECISION. THE ACTUAL TYPE 

C IS AS DECLARED. 

C 

C ERROR CODES . 

C 0 - NO ERROR. 

C 2007 - RECORD TYPE IN THE DATASET IS NOT REAL. 

C 2008 . ERROR IN GMGETN DETECTED BY LMERCD. 

C 

C GAL-PROCESSOR ENTRY POINTS . 

C GMCORN, GMGETN, LMERCD, EMSG. 

C 

Q**«* ******** I*************** **************************************************** 

C 

SUBROUTINE GTRECF ( RECNUM, FBUF, LEN ) 

C 

INTEGERS RECNUM, LEN 
DOUBLE PRECISION FBUF(l) 

1****************************************************************************** 


CHARACTER* 40 DUMMY1 
CHARACTBR*B1 CDUMMY 
CHARACTER* 4 RTYPB 

INTEGER*4 IDSN , LDI , NLEN , NREC , TRACE 
INTEGER*4 MSGLVL, IERR, DUMMY2 

COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE, 

1 TRACE 

COMMON /CSMUSR/ DUMMY1 , MSGLVL, IERR, DUMMY2 , CDUMMY(7) 
INTEGER*4 LMERCD 

:****************************************************************************** 


CHARACTER* 4 BUFTYP 
CHARACTER* 12 OP, RKEY 
CHARACTER*24 RNAMB 

INTEGBR*4 IERROR, IGAP , IHI , ILO , IOFF , MDIM 


DETECT TYPE MISMATCH 


IF ((RTYPE .NE. ’D ’) .AND. (RTYPE .NE. ’S ’)) GO TO 600 


CONSTRUCT NAME 'RKEY.RECNUM :RE CNUM’ FOR AN INDIVIDUAL RECORD 
MAXIMUM LENGTH IS 24 CHARACTERS: 12 FOR RKEY, 5 FOR EACH 
RECNUM REPRESENTING HIGH AND LOW CYCLES. 


RKEY = ’DATA * 

ILO = RECNUM 
IHI = RECNUM 

CALL GMCORN ( RNAME, RKEY, ILO, IHI ) 


OP ARGUMENT FOR GMGETx: ’MAINKEY/ QUALIFIER’ 
MAXIMUM LENGTH IS 11: 4 FOR KEY AND 6 FOR QUALIFIER 


OP = ’READ /LENGTH ’ 
BUFTYP = ’D ’ 

IGAP = 0 
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IOPP = 0 

CALL GMGETN ( OP, LDI, IDSN, RNAMB, BUFTYP, PBUF, LEN, MDIM, 
1 IGAP, IOPP, TRACE ) 


TEST ERROR CONDITION AFTER AN ERROR- SENSITIVE REFERENCE 
TO THE I/O MANAGER 

IERROR = LMERCD ( IERROR ) 

IF { IERROR .KB. 0 ) GO TO 600 
RETURN 
C 

600 CONTINUE 
IERR a 2007 

IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 
C 

600 CONTINUE 
IERR = 2008 

IF ( MSGLVL .GE. 3 ) CALL EMSG 
RETURN 
C 

END 
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0 ******************************************************************************** 
0 *************************** ********* ******************************************** 
C BMSG ... ERROR MESSAGE HANDLINE ROUTINE 

0 *************************** ********* ********* ************** ********************* 
0 ******************************************************************** ************ 
c 

C PURPOSE - THIS ROUTINE IS USED TO HANDLE ERROR MESSAGES IN 
C SYSTEM-CSM WHICH INTERFACES SPARSPAK-A WITH CSM TESTBED 

C DATABASE. 

C 

C PROGRAM SUBROUTINES - 
C EMSGO, EMSG1, DBMSGO 

C 

£l ***************************************************************************** *** 

c 

SUBROUTINE EMSG 
C 

0 ***************************************** ************************************ *** 

c 

CHARACTER 1 * 1 40 LIBNAM 
CHARACTER* 51 CDUMMY 
INTEGERS IPRNTE, IPRNTS, MAXINT 
INTEGERS MSGLVL, IERR , MAXCSM 
REAL RATIOS, RATIOL, TIME 

C 

0 ******************************************************************************** 

C 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 

C 

INTEGERS LEVEL 
C 

WRITE ( IPRNTE, 11 ) 

11 FORMAT (/5X, ’EMSG - SYSTEM-CSM ERROR ... ’ ) 

C 

c 

C DETERMINE THE TYPE OF MODULE THAT CALLED EMSG, 

C AND CALL THE APPROPRIATE ERROR ROUTINE TO PRINT 

C THE ERROR MESSAGE 

C 

c 

IF ( IERR .GT. 2000 ) GO TO 1000 
C 

LEVEL = (IERR - 1000)/10 + 1 
GO TO ( 100, 200, 300 ) , LEVEL 
C 

100 CONTINUE 

C 

C IERR RANGES FROM 1001 TO 1009 

C 

CALL EMSGO 
RETURN 
C 

200 CONTINUE 

C 

C IERR RANGES FROM 1011 TO 1019 

C 

CALL EMSG1 
RETURN 
C 

300 CONTINUE 

C 

C IERR RANGES FROM 1021 TO 1029 

c 

CALL EMSG2 
RETURN 
C 

1000 CONTINUE 

LEVEL = (IERR - 2000)/l0 + 1 
GO TO ( 1100, 1200 ) , LEVEL 
C 

1100 CONTINUE 

C 

C IERR RANGES FROM 2001 TO 2009 

C 

CALL DEMSGO 
RETURN 
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c 

1200 CONTINUE 
RETURN 
C 

END 
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£ *** ***** WWW** ***** **** ***** ******* * * *** ** ********* * ****** ** ******************* * 
0 ** ** *** *» * *** # ******* * ************** ********* ********* ************** *********** 

C EMSGO ERROR MESSAGES FOR ... 

£**** ********* ***** **** *********************** ************** ******************** 
0 *** * * * ******* ***** *** * * * *w * * * * * * **** ** * ** **** ***** * ******** ****** *** *** * ***** ** 

C 

C PURPOSE - THIS ROUTINE IS AN ERROR MESSAGE PRINTING 
C ROUTINE FOR THE MODULE SPACE. 

C 

0********* **** ***** **** ********** **** * ** ** ****** *** ********* **** ************** ** 

C 

SUBROUTINE EMSGO 

**** ********* ***** ***************************************** **********1********** 

CHARACTER*40 LIBNAM 
CHARACTER*51 CDUMMY 
INTEGER** IPRNTE, IPRNTS, MAXINT 

INTEGER** MSGLVL , IERR , MAXCSM 

INTEGER** BUFMAX, MXUSED, MXREQD, STAGE 
REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 

1 TIME 

COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 
COMMON /CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 

****************************************** t*** ********* ************** *********** 

INTEGER** IERROR 
C 

IERROR = IERR - 1000 
GO TO ( 100, 200 ) , IERROR 
C 

100 CONTINUE 

WRITE ( IPRNTE, 11 ) IERR, STAGE, MXREQD 
11 FORMAT (/10X, 35HSPACE • ERROR NUMBER ,17 

1 /10X, 35HINSUFFICIENT STORAGE . , 

1 /10X, 35HTHB LAST STAGE COMPLETED IS , 17 

1 /10X, 35HTO CONTINUE MAXCSM IS AT LEAST , 17 ) 

RETURN 

C 

200 CONTINUE 
RETURN 
C 

END 
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£•#*******] *********************************************************************** 
0 ********** ************ ********** ********************** ********* ********* ******* 
C EMSG1 ERROR MESSAGES FOR ... 

C ***************************************************************************** ** 
0 *********************************************************** ******************** 
c 

C PURPOSE . THIS ROUTINE IS AN ERROR MESSAGE PRINTING 
C ROUTINE FOR MODULES: LIBOPN, GETJDF, GETDOF 

C 

0 ******************************** ********* ************************************** 

c 

SUBROUTINE EMSGI 
C 

q*************************** **************************************************** 

C 

CHARACTER*40 LIBNAM 
CHARACTER*51 CDUMMY 
INTEGERS IPRNTE, IPRNTS, MAXINT 

INTEGERS MSGLVL , IERR , MAXCSM 

REAL RATIOS, RATIOL, TIME 

C 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 

C 

0 ******************************************************************** *********** 

c 

INTEGERS IERR OR 
C 

IF ( IERR .GT. 1012 ) GO TO 250 
C 

C ERROR FOR SUBROUTINE LIBOPN 

C 

IERR OR = IERR - 1010 
GO TO ( 100, 200 ) , IBRROR 
C 

100 CONTINUE 

C IERR = 1011 

WRITE ( IPRNTE, 11 ) IERR 

11 FORMAT (/ 10X, 35HLIBOPN - ERROR NUMBER , 17 

1 /10X, 35HCANNOT OPEN DATASET LIBRARY. ) 

RETURN 

C 

200 CONTINUE 

C IERR = 1012 

WRITE ( IPRNTE, 22 ) IERR 

22 FORMAT (/10X, 35HLIBOPN - ERROR NUMBER , 17 

1 /10X, 35HMAX LOGICAL DEVICE INDEX = 30 

1 / 1 OX, 35HLDI RETURNED EXCEEDS THIS VALUE. ) 

RETURN 

C 

250 CONTINUE 

IF ( IERR .GT. 1014 ) GO TO 450 
C 

c 

C ERROR FOR SUBROUTINE GETJDF 

C 

IBRROR = IERR - 1012 
GO TO ( 300, 400 ) , IBRROR 
C 

300 CONTINUE 

C IERR = 1013 

C 

WRITE ( IPRNTE, 33 ) IERR 

33 FORMAT (/ 10X, 35HGBTJDF - ERROR NUMBER , IT 

1 / 10X, 3 SHIN CORRECT EXECUTION SEQUENCE. ) 

RETURN 

C 

400 CONTINUE 

C IERR = 1014 
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44 

1 

C 

450 


C 

900 


99 

1 

C 


WRITE ( IPRNTB, 44 ) IERR 

FORMAT ( / 10X, 35HGETJDF - ERROR NUMBER , 17 

/10X, 35HDATASET DOES NOT HAVE ALL DATA, ) 
RETURN 

CONTINUE 

IF ( IERR ,EQ. 1019 ) GO TO 900 
RETURN 

CONTINUE 

IERR = 1019 

WRITE ( IPRNTB, 99 ) IERR 

FORMAT ( /10X, 35HGBTDOF . ERROR NUMBER , 17 

/10X, 35HINCORRECT EXECUTION SEQUENCE. ) 
RETURN 

END 
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C ************* ********* ******************* **** ********* ***** **** ** *** **** ***** ** 

q**** ********* ********************************************** ************* ******* 


C EMSG2 ERROR MESSAGES FOR ... 

C *************************************************************** **************** 

Q ************************************************** **** ***** ******************** 

c 

C PURPOSE - THIS ROUTINE IS AN ERROR MESSAGE PRINTING 
C ROUTINE FOR MODULES: GTZERO, GTCOND, GETIJ, FTFORC, 

C GTMOTI, GTNUM5. 

C 


C ************* ****************************************************************** 
C 

SUBROUTINE EMSG2 
C 

C ************* ****************************************************************** 
C 


c 

1 

c 

0 ************* ****************************************************************** 

c 

INTEGERS IERROR 
C 

IERROR = IERR - 1020 

GO TO ( 100, 200, 300, 400, BOO, 600, 700, 800), IERROR 

c 

100 CONTINUE 

C IERR = 1021 

WRITE ( IPRNTE, 11 ) IERR 

11 FORMAT (/ 10X, 35HGTZERO - ERROR NUMBER , 17 

1 /10X, ’INCORRECT EXECUTION SEQUENCE » ) 

RETURN 

C 

200 CONTINUB 

C IERR = 1022 

WRITE ( IPRNTE, 22 ) IERR 

22 FORMAT (/10X, 35HGTCOND - ERROR NUMBER , 17 

1 / 10X, 'INCORRECT EXECUTION SEQUENCE » ) 

RETURN 

C 

300 CONTINUE 

C IERR = 1023 

WRITE ( IPRNTE, 33 ) IERR 

33 FORMAT (/10X, 35HGETIJ . ERROR NUMBER ,17 

1 /10X, 'INCORRECT EXECUTION SEQUENCE • ) 

RETURN 

C 

400 CONTINUE 

C IERR a* 1024 

WRITE ( IPRNTE, 44 ) IERR 

44 FORMAT (/10X, 3BHGTFORC . ERROR NUMBER ,17 

1 /10X, 'INCORRECT EXECUTION SEQUENCE ' ) 

RETURN 

C 

BOO CONTINUE 

C IERR = 102B 

WRITE ( IPRNTE, BS ) IERR 

BB FORMAT (/10X, 36HGTMOTI . ERROR NUMBER , 17 

1 /10X, 'INCORRECT EXECUTION SEQUENCE * ) 

RETURN 


CHARACTER* 40 LIBNAM 
CHARACTER* B1 CDUMMY 
INTEGERS IPRNTE, IPRNTS, MAXINT 

INTEGERS MSGLVL , IERR , MAXCSM 

REAL RATIOS, RATIOL, TIME 

COMMON / CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, 
TIME 

COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 
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c 

600 CONTINUE 

IERR = 1026 

WRITE ( IPRNTE, 66 ) IERR 

56 FORMAT (/10X, 35HGTMOTI - ERROR NUMBER ,17 

1 / 10X, ’UNEXPECTED NONZERO CONSTRAINT VALUE’ ) 

RETURN 

00 CONTINUE 

IERR = 1027 

WRITE ( IPRNTE, 77 ) IERR 

77 FORMAT (/10X, 35HGTMOTI - ERROR NUMBER ,17 

1 /10X, ’ZERO ENTRY FOR A NONZERO CONSTRAINT OCCURS’) 

RETURN 

00 CONTINUE 

IERR = 1028 

WRITE ( IPRNTE, 88 ) IERR 

88 FORMAT (/10X, 35HGTNUM1 - ERROR NUMBER ,17 

1 /10X, ’INCORRECT EXECUTION SEQUENCE ’ ) 

RETURN 

C 

END 
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q********** ****** *********** ********* ** ************** ********* ******* *********** 
Q ****** * * * He ****** ***** *** *************** * **** DC ** * ** * ******** * ***** * ** ********** * 

C DEMSGO ERROR MESSAGES FOR DATASET ACCESSES 

0 *************************** ********* ********* ***** ************* ********* ******* 
Q ***************************************************************** * * * ** ** * * * * * * * 

c 

C PURPOSE - THIS ROUTINE IS AN ERROR MESSAGE PRINTING 
C FOR MODULES ACCESSING DATASETS. 

C 


q************* ********* ***** ********* ********* ************** ********* ********* ** 

C 

SUBROUTINE DEMSGO 
C 

0************* ********* ************ ********************************************* 

C 

CHARACTER* 40 LIBNAM 
CHARACTERS CDUMMY 
INTEGERS IPRNTE, IPRNTS, MAXINT 

INTEGERS MSGLVL , IERR , MAXCSM 

INTEGERS BUFMAX, MXUSED, MXREQD, STAGE 

REAL RATIOS, RATIOL, TIME 

C 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON f CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, CDUMMY(7) 
COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 
C 

0****************** ****************** ******************************************* 

C 

INTEGERS IERR OR 
C 

IF ( IERR .GT. 2004 ) GO TO 150 
C 

C ERROR FROM SUBROUTINE QKINFO 

C . 

IERROR = IERR - 2000 
GO TO ( 100, 200, 300, 400 ) , IERROR 
C 

100 CONTINUE 

C IERR = 2001 

WRITE ( IPRNTE, 11 ) IERR 

11 FORMAT (/10X, 35HQKINFO - ERROR NUMBER , 17 

1 /10X, 35HLMFIND; CANNOT FIND DATASET. ) 

RETURN 

C 

200 CONTINUE 

C IERR = 2002 

WRITE ( IPRNTE, 22 ) IERR 

22 FORMAT (/10X, 35HQKINFO - ERROR NUMBER , 17 

1 /10X, 35HGMGEKA: RECORD DOES NOT EXIST. ) 

RETURN 

C 

300 CONTINUE 

C IERR = 2003 

WRITE ( IPRNTE, 33 ) IERR 

33 FORMAT (/10X, 35HQKINFO - ERROR NUMBER ,17 

1 /10X, 35HGMGECY: RECORD GROUP KEY UNDEFINED. ) 

RETURN 

C 

400 CONTINUE 

C 

C IERR ss 2004 

WRITE ( IPRNTE, 44 ) IERR 

44 FORMAT(/10X, 35HQKINFO - ERROR NUMBER ,17 

1 /10X, 38HGMGEC Y: SEGMENTED RECORD GROUP NOTED. ) 

RETURN 
C 

450 CONTINUE 

C 

C ERROR FROM SUBROUTINE GETRECI OR GTRECF 
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IBRROR = IERR - 2004 

GO TO ( 500, 600, 700, 800, 900 ) IERROR 
C 

600 CONTINUE 
C IERR = 2005 

WRITE ( IPRNTE, 55 ) IERR 

55 FORMAT(/10X, 35HGETRECI - ERROR NUMBER ,17 

1 /10X, 35HRECORD TYPE MISMATCH ... ) 

RETURN 
C 

600 CONTINUE 

C IERR = 2006 

C 

WRITE ( IPRNTE, 66 ) IERR 

66 FORMAT(/10X, 36HGETRECI - ERROR NUMBER ,17 

1 /10X, 35HGMGETN: ERROR DETECTED BY LMERCD... ) 

RETURN 
C 

700 CONTINUE 

C IERR = 2007 

WRITE ( IPRNTE, 77 ) IERR 

77 FORMAT(/10X, 36HGETRECF - ERROR NUMBER ,17 

1 /10X, 35HRECORD TYPE MISMATCH ... ) 

RETURN 
C 

800 CONTINUE 

C IERR = 2008 

C 

WRITE ( IPRNTE, 88 ) IERR 

88 FORMAT(/10X, 35HGETRECF - ERROR NUMBER ,17 

1 /10X, 36HGMGBTN: ERROR DETECTED BY LMERCD... ) 

RETURN 
C 

900 CONTINUE 

C IERR = 2009 

WRITE ( IPRNTE, 99 ) IERR, BUFMAX 
99 FORMAT(/10X, 35HQKINFO . ERROR NUMBER ,17 

1 /10X, 35HBUFMAX MUST BE AT LEAST ,17) 

RETURN 
C 

END 
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* *** ***** **** **** ***** ***** 4c* *** ********* **** * ** ** *** ****** *** * ***** **** ***** ** 
********* * ** * ** * **** 4> *** * * * * * ****** 4c ****** * * * * ************ ****** * * ************ * 

GETSOL RETRIEVE TESTBED SOLUTION ... 

************************************************** ***************************** 
******************************************************************************* 

PURPOSE - RETRIEVE THE TESTBED SOLUTION. ASSUMING THAT THE TESTBED 
SOLUTION IS CORRECT, THE MAXIMUM RELATIVE ERROR IS THEN COMPUTED 
FOR EACH COMPOMENT IN THE SOLUTION VECTOR RETURNED BY SPARSPAK-A 
SOLVER “SOLVES”. 

INPUT PARAMETERS - 

SOL - THE LEADING NEQNS LOCATIONS OF THIS VECTOR CONTAIN 
THE SOLUTION RETURNED BY SPARSPAK-A LINEAR SYSTEM 
SOLVER. 

WORKING PARAMETER . 

FBUF - A REAL OR DOUBLE PRECISION BUFFER OF SIZE BUFMAX. 

THE ACTUAL TYPE IS AS DECLARED. 

OUTPUT PARAMETERS - 

RATIO - THE MAXIMUM RELATIVE ERROR ENCOUNTERED. 
******************************************************************************* 


SUBROUTINE GETSOL ( FBUF, SOL, RATIO ) 

DOUBLE PRECISION FBUF(l), SOL(l), RATIO 
******************************************************************************* 


CHARACTER*^ LIBNAM 

CHARACTER*51 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 
CHARACTER* 4 RTYPE 

INTEGER** IPRNTE, IPRNTS, MAXINT 

INTEGER** IDSN , LD1 , NLEN , NREC .TRACE 
INTEGER** MSGLVL, IBRR , MAXCSM 

INTEGER** MAXDOF , NEQNS , NUMJNT 

REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMSPK/ IDSN , LDI , NLEN , NREC , RTYPE , 

1 TRACE 

COMMON / CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLF, APPLM, STATD 

COMMON /PRBLEM/ MAXDOF , NEQNS , NUMJNT 

******************************************************************************** 


INTEGER** I, II, LEN, NITEMS, INDEX, MAXIND 
DOUBLE PRECISION DELTAX, CSM, WEHAVE, CSMMAX 

WRITE ( IPRNTS, 11 ) 

11 FORMAT (/ 5X, 40HGETSOL - COMPARE WITH TESTBED SOLN ... ) 

IF ( IERR .NB. 0 ) GO TO 300 


ACCESS RECORDS IN DATA SET 'STAT.DISP.* » 
TO RETRIEVE NEQNS SOLUTIONS 


CALL QKINFO ( STATD ) 

IF ( IERR .NE. 0 ) GO TO 999 

TRACE = TRACE + 10 

RATIO = 0.0D0 

NITEMS = 0 

CSMMAX as 0.0D0 

DO 100 I = 1, NREC 

LEN = MIN0 ( NEQNS - NITEMS, NLEN ) 
IF ( LEN .GT. 0 ) THEN 

READ NEXT RECORD 


CALL GTRBCF ( I, FBUF, LEN ) 
IF ( IERR .NB. 0 ) RETURN 


COMPUTE THE MAXIMUM RELATIVE ERROR 


131 


C FBUP CONTAINS THE DATABASE SOLUTION 

C 

DO 200 II = 1, LEN 

NITEMS = NITEMS + 1 
C 

C GET THE COMPONENT WITH MAXIMUM MAGNITUDE 

C 

IF ( DABS (FBUF(II)) .GT. CSMMAX ) THEN 
CSMMAX = DABS (FBUF(II)) 

MAXIND = NITEMS 
ENDIF 

DELTAX = DABS ( FBUF(II) • SOL(NITEMS) ) 

IF ( FBUF(II) .NE. 0.0D0 ) 

1 DELTAX = DELTAX /DAB S(FBUF(II)) 

IF ( DELTAX .GT. RATIO ) THEN 
RATIO = DELTAX 
INDEX = NITEMS 

C 

C SAVE THE PAIR WHICH CAUSES MAX REL ERR 

C 

CSM = FBUF(II) 

WEHAVE = S OL(INDEX) 

ENDIF 

200 CONTINUE 

ENDIF 

100 CONTINUE 

C 

C SUMMARY 

C 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 21 ) STATD, RATIO, 

1 INDEX, CSM, WEHAVE 

21 FORMAT( /I OX, ’MAX. REL ERR COMPARED TO », A51, 

1 /10X, ’IS El 4.7, * IN COMPONENT’, IB, 

1 / 10X, ’CSM SOL = ’, B21.14, » WE HAVE \E21.14 ) 

RETURN 

C 

300 CONTINUE 

C 

C ERROR HANDLING .... (NOT INCLUDED IN EMSG) 

C 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 31 ) 

31 FORMAT (/10X, 3BHGETS OL.INC ORRE CT EXECUTION SEQUENCE ) 
RETURN 
C 

999 CONTINUE 

C 

C ERROR HANDLING .... (NOT INCLUDED IN EMSG) 

C 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 91 ) STATD 
91 FORMAT( / 10X, ’CANNOT FIND AB1 ) 

RETURN 

C 

END 
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**** mm******* ***** **** ***** *«*»* **** ***** **** ***** **** ***** **** ***** **** ***** *** 
************* ********* ********** **** ********* ***** **** ***** **** * ******** ******** 
STATCS PRINT STATISTICS 

******************************************************************************** 
* ******************************************* ****** **** ***** ****** *** **** ******** 

PURPOSE - THIS ROUTINE PRINTS TIME AND STORAGE REQUIREMENTS OP 
THE CURRENT RUN. 

******************************************************************************** 


SUBROUTINE STATCS 

CHARACTER* 40 LIBNAM 

CHAR A CTER*51 JDFSET, KMAP, KSPAR, CON, APPLE, APPLM, STATD 

INTEGER*! IPRNTE, IPRNTS, MAXINT 

INTEGERS MSGLVL , IERR , MAXCSM 

INTEGERM DOF, BUF, MASK, KC, ICLQ, FCON, SPK 

INTEGERS BUFMAX, MXUSED, MXREQD, STAGE 

INTEGER*4 MAXDOF, NEQNS , NUMJNT 

REAL GZTIME, G C TIME , GIJTIM, GFTIMB, GMTIME.GNTIMB, 

1 CSMTIM, CSMSTR 

REAL RATIOS, RATIOL, TIME 

COMMON /CSMSYS/ IPRNTE, IPRNTS, MAXINT, RATIOS, RATIOL, TIME 
COMMON /CSMUSR/ LIBNAM, MSGLVL, IERR , MAXCSM, 

1 JDFSET, KMAP, KSPAR, CON, APPLE, APPLM, STATD 

COMMON / C SMMAP / DOF, BUF, MASK, KC, ICLQ, FCON, SPK 
COMMON / CSMCON/ BUFMAX, MXUSED, MXREQD, STAGE 

COMMON /CSMDTA / GZTIME, GCTIME, GIJTIM, GFTIME, GMTIME.GNTIME, 

1 CSMTIM, CSMSTR 

COMMON /PRBLEM/ MAXDOF, NEQNS , NUMJNT 

******************************************************************************** 
WRITE ( IPRNTS, 11 ) 

11 FORMAT (/5X, 40HSTATCS . SYSTEM-CSM STATISTICS ... ) 

C 

IF ( STAGE .GE. 20 ) GO TO 100 
WRITE (IPRNTS.22) 

22 FORMAT (/10X, 35HNO STATISTICS AVAILABLE. ) 

RETURN 

C 

100 CONTINUE 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 33 ) MAXCSM 
33 FORMAT (/10X, 35HSIZE OF STORAGE ARRAY (MAXCSM) , 110 ) 

IF ( MSGLVL .GE. 2 ) WRITE ( IPRNTS, 44 ) NUMJNT.MAXDOF, NEQNS 

44 FORMAT (/10X, 35HNUMBER OF JOINTS ,110 

1 /10X, 35HMAX DEGREE OF FREEDOME PER JOINT , 110 

1 /10X, 35HNUMBER OF EQUATIONS , 110 ) 

IF ( MSGLVL .GE. 3 ) THEN 
WRITE ( IPRNTS, 45 ) 

WRITE ( IPRNTS, 46 ) D OF, BUF .MASK ,K C,ICLQ ,FC ON, SPK 

45 FORMAT (/10X, 35HADDRESSES OF ARRARYS ) 

46 FORMAT (/10X, 10HDOF ,110 

1 /10X, 10HBUF ,110 

1 /10X, 10HMASK , 110 

1 /10X, 10HKC ,110 

1 /10X, 10HICLQ ,110 

1 /10X, 10HFCON ,110 

1 /10X, 10HSPK , 110 ) 

ENDIF 

c 

CSMSTR = MXREQD 

WRITE (IPRNTS, 133) CSMTIM, CSMSTR 
133 FORMAT ( 10X, 35HTOTAL CSM-TIME REQUIRED , F13.3 

1 / 10X, 35HMAXIMUM CSM-STORAGE REQUIRED , F10.0 ) 

RETURN 

C 

END 
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